[cig-commits] r6288 -
short/3D/PyLith/branches/pylith-0.8/pylith3d/libpylith3d
willic3 at geodynamics.org
willic3 at geodynamics.org
Mon Mar 19 09:56:29 PDT 2007
Author: willic3
Date: 2007-03-19 09:56:29 -0700 (Mon, 19 Mar 2007)
New Revision: 6288
Removed:
short/3D/PyLith/branches/pylith-0.8/pylith3d/libpylith3d/lubksb.f
short/3D/PyLith/branches/pylith-0.8/pylith3d/libpylith3d/ludcmp.f
Modified:
short/3D/PyLith/branches/pylith-0.8/pylith3d/libpylith3d/Makefile.am
short/3D/PyLith/branches/pylith-0.8/pylith3d/libpylith3d/mat_1.f
short/3D/PyLith/branches/pylith-0.8/pylith3d/libpylith3d/mat_5.f
short/3D/PyLith/branches/pylith-0.8/pylith3d/libpylith3d/mat_6.f
short/3D/PyLith/branches/pylith-0.8/pylith3d/libpylith3d/mat_7.f
short/3D/PyLith/branches/pylith-0.8/pylith3d/libpylith3d/mat_8.f
short/3D/PyLith/branches/pylith-0.8/pylith3d/libpylith3d/mat_9.f
short/3D/PyLith/branches/pylith-0.8/pylith3d/libpylith3d/matmod_def.f
short/3D/PyLith/branches/pylith-0.8/pylith3d/libpylith3d/newt.f
Log:
Completely rearranged material models to correspond with current
preferred methods, although more testing is needed.
Also removed Numerical Recipes linear algebra routines and replaced
them with LAPACK calls.
Modified: short/3D/PyLith/branches/pylith-0.8/pylith3d/libpylith3d/Makefile.am
===================================================================
--- short/3D/PyLith/branches/pylith-0.8/pylith3d/libpylith3d/Makefile.am 2007-03-19 16:50:26 UTC (rev 6287)
+++ short/3D/PyLith/branches/pylith-0.8/pylith3d/libpylith3d/Makefile.am 2007-03-19 16:56:29 UTC (rev 6288)
@@ -98,8 +98,6 @@
local.f \
localf.f \
localx.f \
- lubksb.f \
- ludcmp.f \
makemsr.F \
mat_1.f \
mat_2.f \
@@ -273,7 +271,7 @@
includes/nunits_dim.inc \
includes/nvisdat_def.inc \
includes/nvisdat_dim.inc \
- includes/parmat_6.inc \
+ includes/parmat_9.inc \
includes/prestr_matinit_ext.inc \
includes/prscal_dim.inc \
includes/rconsts.inc \
Deleted: short/3D/PyLith/branches/pylith-0.8/pylith3d/libpylith3d/lubksb.f
===================================================================
--- short/3D/PyLith/branches/pylith-0.8/pylith3d/libpylith3d/lubksb.f 2007-03-19 16:50:26 UTC (rev 6287)
+++ short/3D/PyLith/branches/pylith-0.8/pylith3d/libpylith3d/lubksb.f 2007-03-19 16:56:29 UTC (rev 6288)
@@ -1,41 +0,0 @@
- SUBROUTINE lubksb(a,n,np,indx,b)
-c
-c... routine to perform backsubstitution on an LU-decomposed matrix.
-c Adapted from Numerical Recipes.
-c
- include "implicit.inc"
-c
-c... subroutine arguments
-c
- integer n,np
- INTEGER indx(n)
- double precision a(np,np),b(n)
-c
-c... local variables
-c
- INTEGER i,ii,j,ll
- double precision sum
-c
- ii=0
- do 12 i=1,n
- ll=indx(i)
- sum=b(ll)
- b(ll)=b(i)
- if (ii.ne.0)then
- do 11 j=ii,i-1
- sum=sum-a(i,j)*b(j)
-11 continue
- else if (sum.ne.0.0d0) then
- ii=i
- endif
- b(i)=sum
-12 continue
- do 14 i=n,1,-1
- sum=b(i)
- do 13 j=i+1,n
- sum=sum-a(i,j)*b(j)
-13 continue
- b(i)=sum/a(i,i)
-14 continue
- return
- END
Deleted: short/3D/PyLith/branches/pylith-0.8/pylith3d/libpylith3d/ludcmp.f
===================================================================
--- short/3D/PyLith/branches/pylith-0.8/pylith3d/libpylith3d/ludcmp.f 2007-03-19 16:50:26 UTC (rev 6287)
+++ short/3D/PyLith/branches/pylith-0.8/pylith3d/libpylith3d/ludcmp.f 2007-03-19 16:56:29 UTC (rev 6288)
@@ -1,73 +0,0 @@
- SUBROUTINE ludcmp(a,n,np,indx,d)
-c
-c... subroutine to perform an LU decomposition on a matrix.
-c Adapted from Numerical Recipes.
-c
- include "implicit.inc"
-c
-c... parameter definitions
-c
- integer NMAX
- double precision TINY
- PARAMETER (NMAX=500,TINY=1.0d-20)
-c
-c... subroutine arguments
-c
- integer n,np
- INTEGER indx(n)
- double precision d,a(np,np)
-c
-c... local variables
-c
- INTEGER i,imax,j,k
- double precision aamax,dum,sum,vv(NMAX)
- d=1.0d0
- do 12 i=1,n
- aamax=0.0d0
- do 11 j=1,n
- if (abs(a(i,j)).gt.aamax) aamax=abs(a(i,j))
-11 continue
- if (aamax.eq.0.0d0) pause 'singular matrix in ludcmp'
- vv(i)=1.0d0/aamax
-12 continue
- do 19 j=1,n
- do 14 i=1,j-1
- sum=a(i,j)
- do 13 k=1,i-1
- sum=sum-a(i,k)*a(k,j)
-13 continue
- a(i,j)=sum
-14 continue
- aamax=0.0d0
- do 16 i=j,n
- sum=a(i,j)
- do 15 k=1,j-1
- sum=sum-a(i,k)*a(k,j)
-15 continue
- a(i,j)=sum
- dum=vv(i)*abs(sum)
- if (dum.ge.aamax) then
- imax=i
- aamax=dum
- endif
-16 continue
- if (j.ne.imax)then
- do 17 k=1,n
- dum=a(imax,k)
- a(imax,k)=a(j,k)
- a(j,k)=dum
-17 continue
- d=-d
- vv(imax)=vv(j)
- endif
- indx(j)=imax
- if(a(j,j).eq.0.0d0)a(j,j)=TINY
- if(j.ne.n)then
- dum=1.0d0/a(j,j)
- do 18 i=j+1,n
- a(i,j)=a(i,j)*dum
-18 continue
- endif
-19 continue
- return
- END
Modified: short/3D/PyLith/branches/pylith-0.8/pylith3d/libpylith3d/mat_1.f
===================================================================
--- short/3D/PyLith/branches/pylith-0.8/pylith3d/libpylith3d/mat_1.f 2007-03-19 16:50:26 UTC (rev 6287)
+++ short/3D/PyLith/branches/pylith-0.8/pylith3d/libpylith3d/mat_1.f 2007-03-19 16:56:29 UTC (rev 6288)
@@ -32,7 +32,7 @@
c Model number: 1
c Model name: IsotropicLinearElastic
c Number material properties: 3
-c Number state variables: 2
+c Number state variables: 12
c Tangent matrix varies with state: False
c Material properties: Density
c Young's modulus
Modified: short/3D/PyLith/branches/pylith-0.8/pylith3d/libpylith3d/mat_5.f
===================================================================
--- short/3D/PyLith/branches/pylith-0.8/pylith3d/libpylith3d/mat_5.f 2007-03-19 16:50:26 UTC (rev 6287)
+++ short/3D/PyLith/branches/pylith-0.8/pylith3d/libpylith3d/mat_5.f 2007-03-19 16:56:29 UTC (rev 6288)
@@ -66,11 +66,12 @@
c
c... local constants
c
- character labelp(4)*15,modelname*49
+ character labelp(4)*15,modelname*37
data labelp/"Density",
& "Young's modulus",
& "Poisson's ratio",
& "Viscosity"/
+ parameter(modelname="Isotropic Linear Maxwell Viscoelastic")
integer mattype
parameter(mattype=5)
c
@@ -81,7 +82,6 @@
cdebug write(6,*) "Hello from mat_prt_5_f!"
c
ierr=0
- modelname="Isotropic Linear Maxwell Viscoelastic"
c
c... output plot results
c
Modified: short/3D/PyLith/branches/pylith-0.8/pylith3d/libpylith3d/mat_6.f
===================================================================
--- short/3D/PyLith/branches/pylith-0.8/pylith3d/libpylith3d/mat_6.f 2007-03-19 16:50:26 UTC (rev 6287)
+++ short/3D/PyLith/branches/pylith-0.8/pylith3d/libpylith3d/mat_6.f 2007-03-19 16:56:29 UTC (rev 6288)
@@ -30,16 +30,26 @@
c... Program segment to define material model 6:
c
c Model number: 6
-c Model name: IsotropicPowerLawMaxwellViscoelastic
-c Number material properties: 5
-c Number state variables: 18
+c Model name: IsotropicLinearGenMaxwellViscoelastic
+c Number material properties: 9
+c Number state variables: 30
c Tangent matrix varies with state: True
c Material properties: Density
-c Young's modulus
-c Poisson's ratio
-c Power-law exponent
-c Viscosity coefficient
+c Young's Modulus
+c Poisson's Ratio
+c Shear Ratio 1
+c Viscosity 1
+c Shear Ratio 2
+c Viscosity 2
+c Shear Ratio 3
+c Viscosity 3
c
+c... Usage notes:
+c For a simple Maxwell model, set one shear ratio to 1.0 and the
+c other two to 0.0.
+c For a standard linear solid, set two shear ratios to 0.0 and the
+c final one to whatever fraction is appropriate.
+c
subroutine mat_prt_6(prop,nprop,matnum,idout,idsk,kw,kp,
& ierr,errstrng)
c
@@ -47,6 +57,7 @@
c
c Error codes:
c 4: Write error
+c 101: Attempt to use undefined material model
c
include "implicit.inc"
c
@@ -62,22 +73,36 @@
c
c... local constants
c
- character labelp(5)*21,modelname*40
+ character labelp(9)*15,modelname*49
data labelp/"Density",
& "Young's modulus",
& "Poisson's ratio",
- & "Power-law exponent",
- & "Viscosity coefficient"/
- parameter(modelname="Isotropic Power-Law Maxwell Viscoelastic")
+ & "Shear Ratio 1",
+ & "Viscosity 1",
+ & "Shear Ratio 2",
+ & "Viscosity 2",
+ & "Shear Ratio 3",
+ & "Viscosity 3"/
+ parameter(modelname=
+ & "Isotropic Linear Generalized Maxwell Viscoelastic")
integer mattype
parameter(mattype=6)
c
c... local variables
c
integer i
+ double precision sfrac
c
cdebug write(6,*) "Hello from mat_prt_6_f!"
c
+ ierr=0
+ sfrac=prop(4)+prop(6)+prop(8)
+ if(sfrac.lt.0.0d0.or.sfrac.gt.1.0d0) then
+ ierr=116
+ errstrng="mat_prt_6"
+ return
+ end if
+c
c... output plot results
c
if(idsk.eq.ione) then
@@ -107,9 +132,9 @@
return
c
700 format(/,5x,"Material number: ",i7,/,5x,
- & "Material type: ",a40,/,5x,
+ & "Material type: ",a37,/,5x,
& "Number of properties: ",i7,/)
- 710 format(15x,a21,3x,1pe15.8)
+ 710 format(15x,a15,3x,1pe15.8)
c
end
c
@@ -167,17 +192,23 @@
subroutine elas_strs_6(prop,nprop,state,state0,ee,scur,dmat,tmax,
& nstate,nstate0,ierr,errstrng)
c
-c... subroutine to compute stresses for the elastic solution. For this
+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 total strain, and a viscous state variable for each Maxwell
+c element. The minimum Maxwell time is computed, even though this
+c is the elastic solution, as an aid in determining the proper time
+c step size for the next step.
c The current total strain is contained in ee and the computed
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 state(13:18) = dimensionless viscous deviatoric strains for
+c element 1
+c state(19:24) = dimensionless viscous deviatoric strains for
+c element 2
+c state(25:30) = dimensionless viscous deviatoric strains for
+c element 3
c
c The state0 array contains initial stresses.
c
@@ -197,32 +228,108 @@
double precision dmat(nddmat),tmax
character errstrng*(*)
c
+c... local constants
+c
+ double precision diag(6),eng(6)
+ data diag/one,one,one,zero,zero,zero/
+ data eng/one,one,one,half,half,half/
+c
c... local variables
c
- double precision e,pr,anpwr,emhu,rmu
- double precision sdev(nstr),sinv1,steff
+ integer i
+ double precision sfrac,e,pr,vise,rmu,rmue,tmaxe,emean
+ double precision edev(6)
c
cdebug write(6,*) "Hello from elas_strs_6_f!"
c
+ ierr=0
+ sfrac=prop(4)+prop(6)+prop(8)
+ if(sfrac.lt.0.0d0.or.sfrac.gt.1.0d0) then
+ ierr=116
+ errstrng="elas_strs_6"
+ return
+ end if
+c
+c... compute elastic stresses assuming material matrix has already
+c been computed.
+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 strain quantities to initialize viscous strain
+c
+ emean=third*(ee(1)+ee(2)+ee(3))
+ do i=1,nstr
+ edev(i)=ee(i)*eng(i)-emean*diag(i)
+ end do
+ do i=1,3
+ call dcopy(nstr,edev,ione,state(nstr*(i-1)+13),ione)
+ end do
+c
c... compute Maxwell time for current stress state
c
e=prop(2)
pr=prop(3)
- anpwr=prop(4)
- emhu=prop(5)
rmu=half*e/(one+pr)
- call invar(sdev,sinv1,steff,scur)
tmax=big
- if(steff.ne.zero) tmax=((emhu/steff)**(anpwr-one))*emhu/rmu
+ do i=1,3
+ rmue=prop(2*(i-1)+4)*rmu
+ vise=prop(2*(i-1)+5)
+ if(rmue.ne.zero) tmax=min(tmax,vise/rmue)
+ end do
return
end
c
c
+ function compdq_6(deltp,tau)
+c
+c... function to compute the viscous state variable increment delta-q.
+c
+ include "implicit.inc"
+c
+c... parameter definitions
+c
+ include "ndimens.inc"
+ include "rconsts.inc"
+ integer nterms
+ double precision tfrac
+ parameter(nterms=5,tfrac=1.0d-5)
+c
+c... subroutine arguments
+c
+ double precision compdq_6,deltp,tau
+c
+c... intrinsic functions
+c
+ intrinsic exp
+c
+c... local variables
+c
+ integer i
+ double precision fsign,fac,frac
+c
+ compdq_6=zero
+ if(tau.eq.zero) return
+ if(tau.lt.tfrac*deltp) then
+ fsign=one
+ fac=one
+ frac=one
+ compdq_6=one
+ do i=2,nterms
+ fac=fac*dble(i)
+ fsign=-one*fsign
+ frac=frac*(deltp/tau)
+ compdq_6=compdq_6+fsign*frac/fac
+ end do
+ else
+ compdq_6=tau*(one-exp(-deltp/tau))/deltp
+ end if
+ return
+ end
+c
+c
subroutine td_matinit_6(state,dstate,state0,dmat,prop,rtimdat,
& rgiter,iopt,ntimdat,iddmat,tmax,nstate,nstate0,nprop,matchg,ierr,
& errstrng)
@@ -241,7 +348,6 @@
c... parameter definitions
c
include "ndimens.inc"
- include "nconsts.inc"
include "rconsts.inc"
c
c... subroutine arguments
@@ -259,11 +365,16 @@
include "rgiter_dim.inc"
include "ntimdat_dim.inc"
c
+c... external functions
+c
+ double precision compdq_6
+ external compdq_6
+c
c... local variables
c
- double precision e,pr,anpwr,emhu,rmu
- double precision sdev(nstr),sinv1,steff
- double precision cmat(nddmat)
+ integer i
+ double precision sfrac,frac0,e,pr,rmu,rmue,vise,bulkm,sfac,shfac
+ double precision rmut,tau
c
c... included variable definitions
c
@@ -273,223 +384,61 @@
c
cdebug write(6,*) "Hello from td_matinit_6_f!"
c
-c
-c... get stress invariants and compute elastic material matrix
-c
- if(iopt.eq.1) then
- call invar(sdev,sinv1,steff,state)
- else
- call invar(sdev,sinv1,steff,dstate)
- end if
- call elas_mat_6(dmat,prop,iddmat,nprop,ierr,errstrng)
- tmax=big
-c
-c... if second deviatoric invariant is zero, use only elastic solution
-c
- if(steff.eq.zero) then
+ ierr=0
+ sfrac=prop(4)+prop(6)+prop(8)
+ if(sfrac.lt.0.0d0.or.sfrac.gt.1.0d0) then
+ ierr=116
+ errstrng="td_matinit_6"
return
-c
-c... otherwise, invert elastic matrix and augment it by the viscous
-c Jacobian.
-c
- else
-c
-c... invert elastic matrix
-c
- call invsymp(nddmat,nstr,dmat,ierr,errstrng)
- if(ierr.ne.izero) return
-c
-c... define material properties
-c
- e=prop(2)
- pr=prop(3)
- anpwr=prop(4)
- emhu=prop(5)
- rmu=half*e/(one+pr)
-c
-c... compute viscous Jacobian and add it to inverted elastic matrix
-c
- call dbds_6(cmat,nddmat,iddmat,sdev,nstr,steff,anpwr,emhu,
- & alfap,deltp)
- call daxpy(nddmat,one,cmat,ione,dmat,ione)
-c
-c... invert augmented matrix
-c
- call invsymp(nddmat,nstr,dmat,ierr,errstrng)
- if(ierr.ne.izero) return
- tmax=((emhu/steff)**(anpwr-one))*emhu/rmu
end if
- return
- end
c
+c... get material properties
c
- subroutine jaccmp_6(scur,n,fvec,nfjsize,fjac,rpar,nrpar,ipar,
- & nipar,ierr,errstrng)
+ call fill(dmat,zero,nddmat)
+ e=prop(2)
+ pr=prop(3)
+ rmu=half*e/(one+pr)
+ bulkm=third*e/(one-two*pr)
c
-c... subroutine to form the Jacobian used in finding the zero of the
-c stress function.
-c** Note that in this version fjac is always assumed to be a symmetric
-c matrix stored in packed format.
+c... compute Prony series term
c
- include "implicit.inc"
+ frac0=one-sfrac
+ shfac=frac0
+ tmax=big
+ do i=1,3
+ rmue=prop(2*(i-1)+4)
+ rmut=rmu*rmue
+ vise=prop(2*(i-1)+5)
+ tau=zero
+ if(rmue.ne.zero) then
+ tau=vise/rmut
+ tmax=min(tmax,tau)
+ shfac=shfac+rmue*compdq_6(deltp,tau)
+ end if
+ end do
c
-c... parameter definitions
+c... compute stress-strain matrix
c
- include "ndimens.inc"
- include "nconsts.inc"
- include "rconsts.inc"
- include "parmat_6.inc"
-c
-c... subroutine arguments
-c
- integer n,nfjsize,nrpar,nipar,ierr
-c*** note that the dimension below is a bit kludgy and does not allow
-c*** anything else to be put in the ipar array
- integer ipar(nstr,nstr)
- double precision scur(n),fvec(n),fjac(nfjsize),rpar(nrpar)
- character errstrng*(*)
-c
-c... local variables
-c
- integer i,j
-ctest integer iddmat(nstr,nstr)
-ctest equivalence(ipar(indiddmat),iddmat)
- double precision cmat(nddmat)
- double precision sdev(nstr),sinv1,steff
-c
-cdebug write(6,*) "Hello from jaccmp_6_f!"
-c
-c
-c... get stress invariants and a copy of inverted elasticity matrix
-c
- call invar(sdev,sinv1,steff,scur)
- call dcopy(nddmat,rpar(inddmati),ione,fjac,ione)
-c
-c... if second deviatoric invariant is zero, use only elastic solution
-c
- if(steff.eq.zero) then
- return
-c
-c... otherwise, augment inverted elastic matrix by the viscous Jacobian.
-c
- else
-c
-c... compute viscous Jacobian and add it to inverted elastic matrix
-c
- call dbds_6(cmat,nddmat,ipar,sdev,nstr,steff,
- & rpar(indanpwr),rpar(indemhu),rpar(indalfap),rpar(inddeltp))
- call daxpy(nddmat,one,cmat,ione,fjac,ione)
- end if
+ shfac=third*rmu*shfac
+ dmat(iddmat(1,1))=bulkm+four*shfac
+ dmat(iddmat(2,2))=dmat(iddmat(1,1))
+ dmat(iddmat(3,3))=dmat(iddmat(1,1))
+ dmat(iddmat(1,2))=bulkm-two*shfac
+ dmat(iddmat(1,3))=dmat(iddmat(1,2))
+ dmat(iddmat(2,3))=dmat(iddmat(1,2))
+ dmat(iddmat(4,4))=three*shfac
+ dmat(iddmat(5,5))=dmat(iddmat(4,4))
+ dmat(iddmat(6,6))=dmat(iddmat(4,4))
return
end
c
c
- subroutine dbds_6(cmat,nddmat,iddmat,sdev,nstr,steff,anpwr,emhu,
- & alfap,deltp)
-c
-c... subroutine to compute derivatives of viscous strain rate with
-c respect to stress.
-c
- implicit none
-c
-c... parameter definitions
-c
- include "nconsts.inc"
- include "rconsts.inc"
-c
-c... subroutine arguments
-c
- integer nddmat,nstr
- integer iddmat(nstr,nstr)
- double precision cmat(nddmat),sdev(nstr),steff,anpwr,emhu
- double precision alfap,deltp
-c
-c... local variables
-c
- double precision sigma,sigmap,sxx,syy,szz,sxy,syz,sxz
-c
-cdebug write(6,*) "Hello from dbds_6!"
-c
- call fill(cmat,zero,nddmat)
- sigma=(steff/emhu)**(anpwr-one)
- sigmap=(anpwr-one)*sigma
- sigma=half*alfap*deltp*sigma/emhu
- sigmap=fourth*alfap*deltp*sigmap/emhu
- sxx=sdev(1)/steff
- syy=sdev(2)/steff
- szz=sdev(3)/steff
- sxy=two*sdev(4)/steff
- syz=two*sdev(5)/steff
- sxz=two*sdev(6)/steff
-c
-c... compute viscous Jacobian
-c
- cmat(iddmat(1,1))=sigma*two*third+sigmap*sxx*sxx
- cmat(iddmat(1,2))=-sigma*third +sigmap*sxx*syy
- cmat(iddmat(1,3))=-sigma*third +sigmap*sxx*szz
- cmat(iddmat(1,4))= sigmap*sxx*sxy
- cmat(iddmat(1,5))= sigmap*sxx*syz
- cmat(iddmat(1,6))= sigmap*sxx*sxz
- cmat(iddmat(2,2))=sigma*two*third+sigmap*syy*syy
- cmat(iddmat(2,3))=-sigma*third +sigmap*syy*szz
- cmat(iddmat(2,4))= sigmap*syy*sxy
- cmat(iddmat(2,5))= sigmap*syy*syz
- cmat(iddmat(2,6))= sigmap*syy*sxz
- cmat(iddmat(3,3))=sigma*two*third+sigmap*szz*szz
- cmat(iddmat(3,4))= sigmap*szz*sxy
- cmat(iddmat(3,5))= sigmap*szz*syz
- cmat(iddmat(3,6))= sigmap*szz*sxz
- cmat(iddmat(4,4))=sigma*two +sigmap*sxy*sxy
- cmat(iddmat(4,5))= sigmap*sxy*syz
- cmat(iddmat(4,6))= sigmap*sxy*sxz
- cmat(iddmat(5,5))=sigma*two +sigmap*syz*syz
- cmat(iddmat(5,6))= sigmap*syz*sxz
- cmat(iddmat(6,6))=sigma*two +sigmap*sxz*sxz
- return
- end
-c
-c
- subroutine betacmp_6(strs,beta,sdev,seff,sinv1,emhu,anpwr,nstr)
-c
-c... subroutine to compute viscous strain rate vector beta
-c Routine also returns stress invariants
-c
- implicit none
-c
-c... parameter definitions
-c
- include "nconsts.inc"
- include "rconsts.inc"
-c
-c... subroutine arguments
-c
- integer nstr
- double precision strs(nstr),beta(nstr),sdev(nstr),seff,sinv1
- double precision emhu,anpwr
-c
-c... local variables
-c
- double precision rm
-c
-c... compute stress invariants and constants for loop
-c
- call fill(beta,zero,nstr)
- call invar(sdev,sinv1,seff,strs)
- rm=half*(seff/emhu)**(anpwr-one)/emhu
-c
-c... compute strain rate components
-c
- call daxpy(nstr,rm,sdev,ione,beta,ione)
- return
- end
-c
-c
subroutine td_strs_6(state,dstate,state0,ee,scur,dmat,prop,
& rtimdat,rgiter,iter,ntimdat,iddmat,tmax,nstate,nstate0,nprop,
& matchg,ierr,errstrng)
c
-c... subroutine to compute the current values for stress, total strain,
-c and viscous strain increment for the time-dependent solution.
+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
@@ -500,9 +449,6 @@
include "ndimens.inc"
include "nconsts.inc"
include "rconsts.inc"
- include "parmat_6.inc"
- integer nrpar,nipar
- parameter(nrpar=31,nipar=36)
c
c... subroutine arguments
c
@@ -519,24 +465,23 @@
include "rgiter_dim.inc"
include "ntimdat_dim.inc"
c
-c... external functions
+c... external functions
c
- double precision sprod,dasum
- external sprod,funcv_6,jaccmp_6,dasum
+ double precision compdq_6
+ external compdq_6
c
c... local constants
c
- double precision sub
- data sub/-1.0d0/
+ double precision diag(6),eng(6)
+ data diag/one,one,one,zero,zero,zero/
+ data eng/one,one,one,half,half,half/
c
c... local variables
c
- integer i
- double precision e,pr,anpwr,emhu,rmu
- double precision test0,rmt,rmtdt
- double precision sdev(nstr),betat(nstr),betatdt(nstr),sinv1,seff
- double precision rpar(nrpar)
- logical check
+ integer i,j,ind
+ double precision sfrac,e,pr,rmu,bulkm,rmue,rmut,vise,frac0
+ double precision etracetdt,emeantdt,smeantdt,sdevtdt,dq(3)
+ double precision edevtdt,edevt,emeant,deltae,rtime(3)
c
c... included variable definitions
c
@@ -544,137 +489,68 @@
include "rgiter_def.inc"
include "ntimdat_def.inc"
c
-cdebug write(6,*) "Hello from td_strs_6!"
+cdebug write(6,*) "Hello from td_strs_6_f!"
c
ierr=0
- test0=dasum(nstate0,state0,ione)
+ sfrac=prop(4)+prop(6)+prop(8)
+ if(sfrac.lt.zero.or.sfrac.gt.one) then
+ ierr=116
+ errstrng="td_strs_6"
+ return
+ end if
c
-c... define material properties
+c... define material properties and dilatational stress
c
e=prop(2)
pr=prop(3)
- anpwr=prop(4)
- emhu=prop(5)
rmu=half*e/(one+pr)
- tmax=big
- rmt=deltp*(one-alfap)
- rmtdt=deltp*alfap
+ bulkm=third*e/(one-two*pr)
+ etracetdt=ee(1)+ee(2)+ee(3)
+ emeantdt=third*etracetdt
+ smeantdt=bulkm*etracetdt
+ emeant=third*(state(7)+state(8)+state(9))
c
-c... for first iteration, use stresses from previous step as initial
-c guess, otherwise use current estimate.
+c... compute Prony series terms
c
- if(iter.eq.ione) call dcopy(nstr,state,ione,dstate,ione)
+ frac0=one-sfrac
+ tmax=big
+ do i=1,3
+ rmue=prop(2*(i-1)+4)
+ rmut=rmu*rmue
+ vise=prop(2*(i-1)+5)
+ rtime(i)=zero
+ if(rmue.ne.zero) then
+ rtime(i)=vise/rmut
+ tmax=min(tmax,rtime(i))
+ dq(i)=compdq_6(deltp,rtime(i))
+ end if
+ end do
c
-c... compute constant part of iterative solution equation.
+c... compute new stresses and store stress and strain values in dstate
c
-c... current total strain
- call dcopy(nstr,ee,ione,rpar(indconst),ione)
- call dscal(nstr,sub,rpar(indconst),ione)
-c... inverted elasticity matrix times initial stresses
- call elas_mat_6(rpar(inddmati),prop,iddmat,nprop,ierr,errstrng)
- call invsymp(nddmat,nstr,rpar(inddmati),ierr,errstrng)
- if(ierr.ne.izero) return
- if(test0.ne.zero) call dspmv("u",nstr,sub,rpar(inddmati),state0,
- & ione,one,rpar(indconst),ione)
-c... viscous strain from previous time step
- call daxpy(nstr,one,state(13),ione,rpar(indconst),ione)
-c... contribution to viscous strain from stresses at beginning of step
- call betacmp_6(state,betat,sdev,seff,sinv1,emhu,anpwr,nstr)
- call daxpy(nstr,rmt,betat,ione,rpar(indconst),ione)
+ do i=1,nstr
+ edevtdt=eng(i)*ee(i)-diag(i)*emeantdt
+ edevt=eng(i)*state(i+6)-diag(i)*emeant
+ deltae=edevtdt-edevt
+ sdevtdt=frac0*edevtdt
+ do j=1,3
+ rmue=prop(2*(j-1)+4)
+ if(rmue.ne.zero) then
+ ind=i+12+nstr*(j-1)
+ dstate(ind)=state(ind)*exp(-deltp/rtime(j))+dq(j)*deltae
+ sdevtdt=sdevtdt+rmue*dstate(ind)
+ end if
+ end do
+ sdevtdt=two*rmu*sdevtdt
+ dstate(i)=diag(i)*smeantdt+sdevtdt+state0(i)
+ dstate(i+6)=ee(i)
+ scur(i)=dstate(i)
+ end do
c
-c... finish defining parameters for stress computation then call Newton
-c routine to compute stresses.
-c
- rpar(indalfap)=alfap
- rpar(inddeltp)=deltp
- rpar(indemhu)=emhu
- rpar(indanpwr)=anpwr
- call newt(dstate,nstr,rpar,nrpar,iddmat,nipar,funcv_6,jaccmp_6,
- & check,ierr,errstrng)
- if(ierr.ne.izero) return
-c
-c... update state variables for current stress estimates.
-c
- call betacmp_6(dstate,betatdt,sdev,seff,sinv1,emhu,anpwr,nstr)
- call dcopy(nstr,ee,ione,dstate(7),ione)
- call dcopy(nstr,state(13),ione,dstate(13),ione)
- call daxpy(nstr,rmt,betat,ione,dstate(13),ione)
- call daxpy(nstr,rmtdt,betatdt,ione,dstate(13),ione)
- call dcopy(nstr,dstate,ione,scur,ione)
-c
-c... compute current Maxwell time
-c
- call invar(sdev,sinv1,seff,dstate)
- if(seff.ne.zero) tmax=((emhu/seff)**(anpwr-one))*emhu/rmu
-c
return
end
c
c
- subroutine funcv_6(n,scur,r,rpar,nrpar,ipar,nipar,
- & ierr,errstrng)
-c
-c... routine to compute the vector of functions(r) evaluated at the
-c given stress state (contained in scur).
-c
-c... The rpar array contains:
-c 1-6: Constant part of vector, which includes current total
-c strain, viscous strain from previous step, stress-related
-c quantities from prefious time step, and prestress term.
-c 7-27: Inverted elastic material matrix.
-c 28: Integration parameter alfap.
-c 29: Time step size deltp.
-c 30: Viscosity coefficient emhu.
-c 31: Power-law exponent anpwr.
-c
-c... The ipar array contains:
-c 1-36: The iddmat index array
-c
- implicit none
-c
-c... parameter definitions
-c
- include "ndimens.inc"
- include "nconsts.inc"
- include "rconsts.inc"
- include "parmat_6.inc"
-c
-c... subroutine arguments
-c
- integer n,nrpar,nipar,ierr
- integer ipar(nipar)
- double precision scur(n),r(n),rpar(nrpar)
- character errstrng*(*)
-c
-c... local data
-c
- double precision sub
- data sub/-1.0d0/
-c
-c... local variables
-c
- double precision betatdt(nstr),rmtdt
- double precision sdev(nstr),seff,sinv1
-c
-c... viscous strain rate multiplier
-c
- rmtdt=rpar(indalfap)*rpar(inddeltp)
-c
-c... compute current viscous strain rate
-c
- call betacmp_6(scur,betatdt,sdev,seff,sinv1,
- & rpar(indemhu),rpar(indanpwr),nstr)
-c
-c... compute contributions to r and accumulate
-c
- call dcopy(nstr,rpar(indconst),ione,r,ione)
- call daxpy(nstr,rmtdt,betatdt,ione,r,ione)
- call dspmv("u",nstr,one,rpar(inddmati),scur,ione,one,r,ione)
-c
- return
- end
-c
-c
subroutine td_strs_mat_6(state,dstate,state0,ee,scur,dmat,prop,
& rtimdat,rgiter,iter,ntimdat,iddmat,tmax,nstate,nstate0,nprop,
& matchg,ierr,errstrng)
@@ -707,15 +583,9 @@
include "rgiter_dim.inc"
include "ntimdat_dim.inc"
c
-c... external routines
-c
-c
-c... local constants
-c
-c
c... local variables
c
- integer iopt
+ integer iopt
c
c... included variable definitions
c
@@ -723,19 +593,15 @@
include "rgiter_def.inc"
include "ntimdat_def.inc"
c
- ierr=0
- iopt=2
+cdebug write(6,*) "Hello from td_strs_mat_6_f!
c
-c... compute current stress estimates for time-dependent solution
-c
- call td_strs_6(state,dstate,state0,ee,scur,dmat,prop,
- & rtimdat,rgiter,iter,ntimdat,iddmat,tmax,nstate,nstate0,nprop,
- & matchg,ierr,errstrng)
-c
-c... compute tangent stiffness corresponding to new stress estimates
-c
- call td_matinit_6(state,dstate,state0,dmat,prop,rtimdat,
- & rgiter,iopt,ntimdat,iddmat,tmax,nstate,nstate0,nprop,matchg,ierr,
+ iopt=2
+ call td_matinit_6(state,dstate,state0,dmat,prop,rtimdat,rgiter,
+ & iopt,ntimdat,iddmat,tmax,nstate,nstate0,nprop,matchg,
+ & ierr,errstrng)
+ if(ierr.ne.izero) return
+ call td_strs_6(state,dstate,state0,ee,scur,dmat,prop,rtimdat,
+ & rgiter,iter,ntimdat,iddmat,tmax,nstate,nstate0,nprop,matchg,ierr,
& errstrng)
c
return
@@ -768,8 +634,10 @@
c
c... local variables
c
- double precision ptmp(10)
+ double precision ptmp(20)
c
+cdebug write(6,*) "Hello from prestr_mat_6_f!"
+c
call dcopy(nprop,prop,ione,ptmp,ione)
if(ipauto.eq.ione) then
ptmp(2)=tyoungs
@@ -800,7 +668,6 @@
c
call dcopy(nstate,state,ione,sout,ione)
call dcopy(nstate,dstate,ione,sout(nstatesmax+ione),ione)
-c
return
end
c
@@ -831,7 +698,7 @@
double precision sub
data sub/-1.0d0/
c
-cdebug write(6,*) "Hello from update_state_6_f!"
+c... local variables
c
call daxpy(nstate,sub,state,ione,dstate,ione)
call daxpy(nstate,one,dstate,ione,state,ione)
Modified: short/3D/PyLith/branches/pylith-0.8/pylith3d/libpylith3d/mat_7.f
===================================================================
--- short/3D/PyLith/branches/pylith-0.8/pylith3d/libpylith3d/mat_7.f 2007-03-19 16:50:26 UTC (rev 6287)
+++ short/3D/PyLith/branches/pylith-0.8/pylith3d/libpylith3d/mat_7.f 2007-03-19 16:56:29 UTC (rev 6288)
@@ -30,26 +30,16 @@
c... Program segment to define material model 7:
c
c Model number: 7
-c Model name: IsotropicLinearGenMaxwellViscoelastic
-c Number material properties: 9
-c Number state variables: 30
+c Model name: IsotropicPowerLawMaxwellViscoelastic
+c Number material properties: 5
+c Number state variables: 18
c Tangent matrix varies with state: True
c Material properties: Density
-c Young's Modulus
-c Poisson's Ratio
-c Shear Ratio 1
-c Viscosity 1
-c Shear Ratio 2
-c Viscosity 2
-c Shear Ratio 3
-c Viscosity 3
+c Young's modulus
+c Poisson's ratio
+c Power-law exponent
+c Viscosity coefficient
c
-c... Usage notes:
-c For a simple Maxwell model, set one shear ratio to 1.0 and the
-c other two to 0.0.
-c For a standard linear solid, set two shear ratios to 0.0 and the
-c final one to whatever fraction is appropriate.
-c
subroutine mat_prt_7(prop,nprop,matnum,idout,idsk,kw,kp,
& ierr,errstrng)
c
@@ -57,7 +47,6 @@
c
c Error codes:
c 4: Write error
-c 101: Attempt to use undefined material model
c
include "implicit.inc"
c
@@ -73,35 +62,22 @@
c
c... local constants
c
- character labelp(9)*15,modelname*49
+ character labelp(5)*21,modelname*40
data labelp/"Density",
& "Young's modulus",
& "Poisson's ratio",
- & "Shear Ratio 1",
- & "Viscosity 1",
- & "Shear Ratio 2",
- & "Viscosity 2",
- & "Shear Ratio 3",
- & "Viscosity 3"/
+ & "Power-law exponent",
+ & "Viscosity coefficient"/
+ parameter(modelname="Isotropic Power-Law Maxwell Viscoelastic")
integer mattype
parameter(mattype=7)
c
c... local variables
c
integer i
- double precision sfrac
c
cdebug write(6,*) "Hello from mat_prt_7_f!"
c
- ierr=0
- modelname="Isotropic Linear Generalized Maxwell Viscoelastic"
- sfrac=prop(4)+prop(6)+prop(8)
- if(sfrac.lt.0.0d0.or.sfrac.gt.1.0d0) then
- ierr=116
- errstrng="mat_prt_7"
- return
- end if
-c
c... output plot results
c
if(idsk.eq.ione) then
@@ -131,9 +107,9 @@
return
c
700 format(/,5x,"Material number: ",i7,/,5x,
- & "Material type: ",a37,/,5x,
+ & "Material type: ",a40,/,5x,
& "Number of properties: ",i7,/)
- 710 format(15x,a15,3x,1pe15.8)
+ 710 format(15x,a21,3x,1pe15.8)
c
end
c
@@ -191,23 +167,17 @@
subroutine elas_strs_7(prop,nprop,state,state0,ee,scur,dmat,tmax,
& nstate,nstate0,ierr,errstrng)
c
-c... subroutine to compute stresses for the elastic solution. For this
+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 a viscous state variable for each Maxwell
-c element. The minimum Maxwell time is computed, even though this
-c is the elastic solution, as an aid in determining the proper time
-c step size for the next step.
+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 stress should be copied to scur.
c
c state(1:6) = Cauchy stress
c state(7:12) = linear strain
-c state(13:18) = dimensionless viscous deviatoric strains for
-c element 1
-c state(19:24) = dimensionless viscous deviatoric strains for
-c element 2
-c state(25:30) = dimensionless viscous deviatoric strains for
-c element 3
+c state(13:18) = viscous strain
c
c The state0 array contains initial stresses.
c
@@ -227,108 +197,32 @@
double precision dmat(nddmat),tmax
character errstrng*(*)
c
-c... local constants
-c
- double precision diag(6),eng(6)
- data diag/one,one,one,zero,zero,zero/
- data eng/one,one,one,half,half,half/
-c
c... local variables
c
- integer i
- double precision sfrac,e,pr,vise,rmu,rmue,tmaxe,emean
- double precision edev(6)
+ double precision e,pr,anpwr,emhu,rmu
+ double precision sdev(nstr),sinv1,steff
c
cdebug write(6,*) "Hello from elas_strs_7_f!"
c
- ierr=0
- sfrac=prop(4)+prop(6)+prop(8)
- if(sfrac.lt.0.0d0.or.sfrac.gt.1.0d0) then
- ierr=116
- errstrng="elas_strs_7"
- return
- end if
-c
-c... compute elastic stresses assuming material matrix has already
-c been computed.
-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 strain quantities to initialize viscous strain
-c
- emean=third*(ee(1)+ee(2)+ee(3))
- do i=1,nstr
- edev(i)=ee(i)*eng(i)-emean*diag(i)
- end do
- do i=1,3
- call dcopy(nstr,edev,ione,state(nstr*(i-1)+13),ione)
- end do
-c
c... compute Maxwell time for current stress state
c
e=prop(2)
pr=prop(3)
+ anpwr=prop(4)
+ emhu=prop(5)
rmu=half*e/(one+pr)
+ call invar(sdev,sinv1,steff,scur)
tmax=big
- do i=1,3
- rmue=prop(2*(i-1)+4)*rmu
- vise=prop(2*(i-1)+5)
- if(rmue.ne.zero) tmax=min(tmax,vise/rmue)
- end do
+ if(steff.ne.zero) tmax=((emhu/steff)**(anpwr-one))*emhu/rmu
return
end
c
c
- function compdq_7(deltp,tau)
-c
-c... function to compute the viscous state variable increment delta-q.
-c
- include "implicit.inc"
-c
-c... parameter definitions
-c
- include "ndimens.inc"
- include "rconsts.inc"
- integer nterms
- double precision tfrac
- parameter(nterms=5,tfrac=1.0d-5)
-c
-c... subroutine arguments
-c
- double precision compdq_7,deltp,tau
-c
-c... intrinsic functions
-c
- intrinsic exp
-c
-c... local variables
-c
- integer i
- double precision fsign,fac,frac
-c
- compdq_7=zero
- if(tau.eq.zero) return
- if(tau.lt.tfrac*deltp) then
- fsign=one
- fac=one
- frac=one
- compdq_7=one
- do i=2,nterms
- fac=fac*dble(i)
- fsign=-one*fsign
- frac=frac*(deltp/tau)
- compdq_7=compdq_7+fsign*frac/fac
- end do
- else
- compdq_7=tau*(one-exp(-deltp/tau))/deltp
- end if
- return
- end
-c
-c
subroutine td_matinit_7(state,dstate,state0,dmat,prop,rtimdat,
& rgiter,iopt,ntimdat,iddmat,tmax,nstate,nstate0,nprop,matchg,ierr,
& errstrng)
@@ -364,16 +258,10 @@
include "rgiter_dim.inc"
include "ntimdat_dim.inc"
c
-c... external functions
-c
- double precision compdq_7
- external compdq_7
-c
c... local variables
c
- integer i
- double precision sfrac,frac0,e,pr,rmu,rmue,vise,bulkm,sfac,shfac
- double precision rmut,tau
+ double precision e,pr,anpwr,emhu,f1,f2,gam,ae,rmu
+ double precision sdev(nstr),sinv1,steff
c
c... included variable definitions
c
@@ -383,49 +271,26 @@
c
cdebug write(6,*) "Hello from td_matinit_7_f!"
c
- ierr=0
- sfrac=prop(4)+prop(6)+prop(8)
- if(sfrac.lt.0.0d0.or.sfrac.gt.1.0d0) then
- ierr=116
- errstrng="td_matinit_7"
- return
- end if
-c
-c... get material properties
-c
call fill(dmat,zero,nddmat)
e=prop(2)
pr=prop(3)
+ anpwr=prop(4)
+ emhu=prop(5)
+ ae=(one+pr)/e
rmu=half*e/(one+pr)
- bulkm=third*e/(one-two*pr)
-c
-c... compute Prony series term
-c
- frac0=one-sfrac
- shfac=frac0
+ f1=third*e/(one-two*pr)
+ call invar(sdev,sinv1,steff,state)
+ gam=half*((steff/emhu)**(anpwr-one))/emhu
tmax=big
- do i=1,3
- rmue=prop(2*(i-1)+4)
- rmut=rmu*rmue
- vise=prop(2*(i-1)+5)
- tau=zero
- if(rmue.ne.zero) then
- tau=vise/rmut
- tmax=min(tmax,tau)
- shfac=shfac+rmue*compdq_7(deltp,tau)
- end if
- end do
-c
-c... compute stress-strain matrix
-c
- shfac=third*rmu*shfac
- dmat(iddmat(1,1))=bulkm+four*shfac
+ if(steff.ne.zero) tmax=((emhu/steff)**(anpwr-one))*emhu/rmu
+ f2=third/(ae+deltp*gam)
+ 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))=bulkm-two*shfac
+ 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))=three*shfac
+ dmat(iddmat(4,4))=half*three*f2
dmat(iddmat(5,5))=dmat(iddmat(4,4))
dmat(iddmat(6,6))=dmat(iddmat(4,4))
return
@@ -436,8 +301,8 @@
& rtimdat,rgiter,iter,ntimdat,iddmat,tmax,nstate,nstate0,nprop,
& matchg,ierr,errstrng)
c
-c... subroutine to compute the current stress for the time-dependent
-c solution.
+c... subroutine to compute the current values for stress, total strain,
+c and viscous strain increment 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
@@ -448,6 +313,8 @@
include "ndimens.inc"
include "nconsts.inc"
include "rconsts.inc"
+ integer nrpar,nipar
+ parameter(nrpar=9,nipar=1)
c
c... subroutine arguments
c
@@ -464,23 +331,29 @@
include "rgiter_dim.inc"
include "ntimdat_dim.inc"
c
-c... external functions
+c... external functions
c
- double precision compdq_7
- external compdq_7
+ double precision sprod
+ external sprod
c
c... local constants
c
- double precision diag(6),eng(6)
+ 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,j,ind
- double precision sfrac,e,pr,rmu,bulkm,rmue,rmut,vise,frac0
- double precision etracetdt,emeantdt,smeantdt,sdevtdt,dq(3)
- double precision edevtdt,edevt,emeant,deltae,rtime(3)
+ integer i
+ integer ipar(nipar)
+ double precision rpar(nrpar)
+ double precision e,pr,anpwr,emhu,rmu,ae,tf
+ double precision emeantdt,smeant,smeantdt,smean0
+ double precision sinv1t,sefft,sinv10,seff0,sefftdt
+ double precision sefftau,gamtau,f1,f2,sdevtdt,sdevtau
+ double precision epptdt(6),sdevt(6),sdev0(6)
+
c
c... included variable definitions
c
@@ -488,68 +361,206 @@
include "rgiter_def.inc"
include "ntimdat_def.inc"
c
-cdebug write(6,*) "Hello from td_strs_7_f!"
+cdebug write(6,*) "Hello from td_strs_7!"
c
ierr=0
- sfrac=prop(4)+prop(6)+prop(8)
- if(sfrac.lt.zero.or.sfrac.gt.one) then
- ierr=116
- errstrng="td_strs_7"
- return
- end if
c
-c... define material properties and dilatational stress
+c... define material properties
c
e=prop(2)
pr=prop(3)
+ anpwr=prop(4)
+ emhu=prop(5)
rmu=half*e/(one+pr)
- bulkm=third*e/(one-two*pr)
- etracetdt=ee(1)+ee(2)+ee(3)
- emeantdt=third*etracetdt
- smeantdt=bulkm*etracetdt
- emeant=third*(state(7)+state(8)+state(9))
+ ae=(one+pr)/e
+ tf=deltp*(one-alfap)
c
-c... compute Prony series terms
+c... define stress and strain quantities
c
- frac0=one-sfrac
- tmax=big
- do i=1,3
- rmue=prop(2*(i-1)+4)
- rmut=rmu*rmue
- vise=prop(2*(i-1)+5)
- rtime(i)=zero
- if(rmue.ne.zero) then
- rtime(i)=vise/rmut
- tmax=min(tmax,rtime(i))
- dq(i)=compdq_7(deltp,rtime(i))
- end if
+ emeantdt=third*(ee(1)+ee(2)+ee(3))
+ smeantdt=e*emeantdt/(one-two*pr)
+ call invar(sdevt,sinv1t,sefft,state)
+ smeant=third*sinv1t
+ call invar(sdev0,sinv10,seff0,state0)
+ smean0=third*sinv10
+ do i=1,nstr
+ epptdt(i)=eng(i)*ee(i)-diag(i)*emeantdt-state(i+12)
end do
c
+c... define parameters needed by effective stress function
+c
+ rpar(1)=ae
+ rpar(2)=anpwr
+ rpar(3)=emhu
+ rpar(4)=sprod(epptdt,epptdt)
+ rpar(5)=two*ae*sprod(epptdt,sdev0)+ae*ae*seff0*seff0
+ rpar(6)=two*tf*(sprod(epptdt,sdevt)+ae*sprod(sdevt,sdev0))
+ rpar(7)=sefft
+ rpar(8)=deltp
+ rpar(9)=alfap
+c
+c... call effective stress driver routine
+c
+ call esfcomp_7(sefftdt,stol,rpar,nrpar,ipar,nipar,ierr,errstrng)
+ if(ierr.ne.izero) return
+c
+c... compute quantities at intermediate time tau used to compute
+c values at end of time step.
+c
+ sefftau=(one-alfap)*sefft+alfap*sefftdt
+ gamtau=half*((sefftau/emhu)**(anpwr-one))/emhu
+ f1=one/(ae+alfap*deltp*gamtau)
+ f2=tf*gamtau
+ tmax=big
+ if(sefftdt.ne.zero) tmax=((emhu/sefftdt)**(anpwr-one))*emhu/rmu
+c
c... compute new stresses and store stress and strain values in dstate
c
do i=1,nstr
- edevtdt=eng(i)*ee(i)-diag(i)*emeantdt
- edevt=eng(i)*state(i+6)-diag(i)*emeant
- deltae=edevtdt-edevt
- sdevtdt=frac0*edevtdt
- do j=1,3
- rmue=prop(2*(j-1)+4)
- if(rmue.ne.zero) then
- ind=i+12+nstr*(j-1)
- dstate(ind)=state(ind)*exp(-deltp/rtime(j))+dq(j)*deltae
- sdevtdt=sdevtdt+rmue*dstate(ind)
- end if
- end do
- sdevtdt=two*rmu*sdevtdt
- dstate(i)=diag(i)*smeantdt+sdevtdt+state0(i)
+ sdevtdt=f1*(epptdt(i)-f2*sdevt(i)+ae*sdev0(i))
+ dstate(i)=sdevtdt+diag(i)*(smeantdt+smean0)
+ scur(i)=dstate(i)
+ sdevtau=(one-alfap)*sdevt(i)+alfap*sdevtdt
dstate(i+6)=ee(i)
- scur(i)=dstate(i)
+ dstate(i+12)=deltp*gamtau*sdevtau
end do
c
return
end
c
c
+ subroutine esfcomp_7(sefftdt,stol,rpar,nrpar,ipar,nipar,
+ & ierr,errstrng)
+c
+c... driver routine to compute the effective stress at the current
+c time step.
+c
+ include "implicit.inc"
+c
+c... parameter definitions
+c
+ include "rconsts.inc"
+c
+c... subroutine arguments
+c
+ integer nrpar,nipar,ierr
+ integer ipar(nipar)
+ double precision sefftdt,stol,rpar(nrpar)
+ character errstrng*(*)
+c
+c... included dimension and type statements
+c
+c
+c... intrinsic functions
+c
+ intrinsic sqrt,max
+c
+c... external routines
+c
+ external esf_7
+ double precision rtsafe
+ external rtsafe
+c
+c... local variables
+c
+ double precision sefft,sefflo,seffhi,xacc
+c
+c... included variable definitions
+c
+c
+cdebug write(6,*) "Hello from esfcomp_7!"
+c
+ ierr=0
+c
+c... use stress from previous step as initial high value, defaulting to
+c stol if this is zero.
+c
+ sefft=rpar(7)
+ seffhi=max(sefft,stol)
+ sefflo=half*seffhi
+ xacc=max(stol*half*(seffhi-sefflo),stol)
+c
+c... bracket the root
+c
+ call zbrac(esf_7,sefflo,seffhi,rpar,nrpar,ipar,nipar,
+ & ierr,errstrng)
+ if(ierr.ne.0) return
+c
+c... compute effective stress
+c
+ sefftdt=rtsafe(esf_7,sefflo,seffhi,xacc,rpar,nrpar,ipar,nipar,
+ & ierr,errstrng)
+ return
+ end
+c
+c
+ subroutine esf_7(seffcur,f,df,rpar,nrpar,ipar,nipar)
+c
+c... subroutine to compute the effective stress function and it's
+c derivative for a given effective stress for a Maxwell power-law
+c viscoelastic material.
+c
+ include "implicit.inc"
+c
+c... parameter definitions
+c
+ include "rconsts.inc"
+c
+c... subroutine arguments
+c
+ integer nrpar,nipar
+ integer ipar(nipar)
+ double precision seffcur,f,df,rpar(nrpar)
+c
+c... included dimension and type statements
+c
+c
+c... local variables
+c
+ double precision ae,anpwr,emhu,a,b,c,sefft,deltp,alfap,d
+ double precision f1,tf
+ double precision gamtau,sefftau,dgam
+c
+c... included variable definitions
+c
+c
+cdebug write(6,*) "Hello from esf_7!"
+c
+c
+c... values from parameter list and a few other constants
+c
+ ae=rpar(1)
+ anpwr=rpar(2)
+ emhu=rpar(3)
+ b=rpar(4)+rpar(5)
+ c=rpar(6)
+ sefft=rpar(7)
+ deltp=rpar(8)
+ alfap=rpar(9)
+ f1=one-alfap
+ tf=deltp*f1
+ d=tf*sefft
+c
+c... additional values dependent on current effective stress
+c
+ sefftau=f1*sefft+alfap*seffcur
+ gamtau=half*((sefftau/emhu)**(anpwr-one))/emhu
+ a=ae+alfap*deltp*gamtau
+c
+c... effective stress function and it's derivative
+c
+ f=a*a*seffcur*seffcur-b+c*gamtau-d*d*gamtau*gamtau
+ dgam=half*alfap*(anpwr-one)*((sefftau/emhu)**(anpwr-two))/
+ & (emhu*emhu)
+ df=two*a*a*seffcur+(two*a*alfap*deltp*seffcur*seffcur+
+ & c-two*d*d*gamtau)*dgam
+cdebug write(6,*) seffcur,f,df
+c
+ return
+ end
+
+c
+c
subroutine td_strs_mat_7(state,dstate,state0,ee,scur,dmat,prop,
& rtimdat,rgiter,iter,ntimdat,iddmat,tmax,nstate,nstate0,nprop,
& matchg,ierr,errstrng)
@@ -566,6 +577,8 @@
include "ndimens.inc"
include "nconsts.inc"
include "rconsts.inc"
+ integer nrpar,nipar
+ parameter(nrpar=9,nipar=1)
c
c... subroutine arguments
c
@@ -582,9 +595,29 @@
include "rgiter_dim.inc"
include "ntimdat_dim.inc"
c
+c... external routines
+c
+ double precision sprod
+ external sprod
+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 iopt
+ integer i
+ integer ipar(nipar)
+ double precision rpar(nrpar)
+ double precision e,pr,anpwr,emhu,rmu,ae,c1,c2,c3,a,c,d,rk1,rk2
+ double precision emeantdt,smeant,smeantdt,smean0
+ double precision sinv1t,sefft,sinv10,seff0,sefftdt
+ double precision sefftau,gamtau,sdevtdt,sdevtau
+ double precision f1,f2,f3,f4,rkdnm,rkfac,dgde,dsde
+ double precision epptdt(6),sdevt(6),sdev0(6)
c
c... included variable definitions
c
@@ -592,17 +625,123 @@
include "rgiter_def.inc"
include "ntimdat_def.inc"
c
-cdebug write(6,*) "Hello from td_strs_mat_7_f!
+c... rather than call td_strs_7, we duplicate the functionality here,
+c since otherwise there would be an enormous amount of duplicate
+c computations.
c
- iopt=2
- call td_matinit_7(state,dstate,state0,dmat,prop,rtimdat,rgiter,
- & iopt,ntimdat,iddmat,tmax,nstate,nstate0,nprop,matchg,
- & ierr,errstrng)
+ ierr=0
+c
+c... define material properties
+c
+ e=prop(2)
+ pr=prop(3)
+ anpwr=prop(4)
+ emhu=prop(5)
+ rmu=half*e/(one+pr)
+ ae=(one+pr)/e
+ c1=one-alfap
+ c2=deltp*c1
+ c3=third*e/(one-two*pr)
+c
+c... define stress and strain quantities
+c
+ emeantdt=ee(1)+ee(2)+ee(3)
+ smeantdt=c3*emeantdt
+ emeantdt=third*emeantdt
+ call invar(sdevt,sinv1t,sefft,state)
+ smeant=third*sinv1t
+ call invar(sdev0,sinv10,seff0,state0)
+ smean0=third*sinv10
+ do i=1,nstr
+ epptdt(i)=eng(i)*ee(i)-diag(i)*emeantdt-state(i+12)
+ end do
+c
+c... define parameters needed by effective stress function
+c
+ rpar(1)=ae
+ rpar(2)=anpwr
+ rpar(3)=emhu
+ rpar(4)=sprod(epptdt,epptdt)
+ rpar(5)=two*ae*sprod(epptdt,sdev0)+ae*ae*seff0*seff0
+ rpar(6)=two*c2*(sprod(epptdt,sdevt)+ae*sprod(sdevt,sdev0))
+ rpar(7)=sefft
+ rpar(8)=deltp
+ rpar(9)=alfap
+c
+c... call effective stress driver routine
+c
+ call esfcomp_7(sefftdt,stol,rpar,nrpar,ipar,nipar,ierr,errstrng)
if(ierr.ne.izero) return
- call td_strs_7(state,dstate,state0,ee,scur,dmat,prop,rtimdat,
- & rgiter,iter,ntimdat,iddmat,tmax,nstate,nstate0,nprop,matchg,ierr,
- & errstrng)
c
+c... compute quantities at intermediate time tau used to compute
+c values at end of time step.
+c
+ sefftau=(one-alfap)*sefft+alfap*sefftdt
+ gamtau=half*((sefftau/emhu)**(anpwr-one))/emhu
+ a=ae+alfap*deltp*gamtau
+ f1=one/a
+ f2=c2*gamtau
+c
+c... compute new stresses and store stress and strain values in dstate
+c
+ do i=1,nstr
+ sdevtdt=f1*(epptdt(i)-f2*sdevt(i)+ae*sdev0(i))
+ dstate(i)=sdevtdt+diag(i)*(smeantdt+smean0)
+ scur(i)=dstate(i)
+ sdevtau=(one-alfap)*sdevt(i)+alfap*sdevtdt
+ dstate(i+6)=ee(i)
+ dstate(i+12)=deltp*gamtau*sdevtau
+ end do
+ tmax=big
+ if(sefftdt.ne.zero) tmax=((emhu/sefftdt)**(anpwr-one))*emhu/rmu
+c
+c... define some constants for tangent computation
+c
+ rk1=half*(anpwr-one)*((sefftau/emhu)**(anpwr-two))/(emhu*emhu)
+ rk2=sefftau-f1*sefft
+ c=rpar(6)
+ d=c2*rpar(7)
+ f3=alfap*deltp*f1
+ f4=gamtau*c2
+ rkdnm=two*a*rk2*(alfap*deltp*rk1*rk2+a)/(alfap*alfap)+
+ & rk1*(c-two*d*d*gamtau)
+ rkfac=rk1/rkdnm
+c
+c... compute tangent matrix
+c
+ dgde=rkfac*(half*epptdt(1)+ae*sdev0(1)-f4*sdevt(1))
+ dsde=f1*(one-dgde*(c2*sdevt(1)+
+ & f3*(epptdt(1)-f4*sdevt(1)+ae*sdev0(1))))
+ dmat(iddmat(1,1))=c3+two*third*dsde
+ dmat(iddmat(1,2))=c3-third*dsde
+ dmat(iddmat(1,3))=dmat(iddmat(1,2))
+c
+ dgde=rkfac*(half*epptdt(2)+ae*sdev0(2)-f4*sdevt(2))
+ dsde=f1*(one-dgde*(c2*sdevt(2)+
+ & f3*(epptdt(2)-f4*sdevt(2)+ae*sdev0(2))))
+ dmat(iddmat(2,2))=c3+two*third*dsde
+ dmat(iddmat(2,3))=c3-third*dsde
+c
+ dgde=rkfac*(half*epptdt(3)+ae*sdev0(3)-f4*sdevt(3))
+ dsde=f1*(one-dgde*(c2*sdevt(3)+
+ & f3*(epptdt(3)-f4*sdevt(3)+ae*sdev0(3))))
+ dmat(iddmat(3,3))=c3+two*third*dsde
+c
+ dgde=rkfac*(half*epptdt(4)+ae*sdev0(4)-f4*sdevt(4))
+ dsde=f1*(one-dgde*(c2*sdevt(4)+
+ & f3*(epptdt(4)-f4*sdevt(4)+ae*sdev0(4))))
+ dmat(iddmat(4,4))=half*dsde
+c
+ dgde=rkfac*(half*epptdt(5)+ae*sdev0(5)-f4*sdevt(5))
+ dsde=f1*(one-dgde*(c2*sdevt(5)+
+ & f3*(epptdt(5)-f4*sdevt(5)+ae*sdev0(5))))
+ dmat(iddmat(5,5))=half*dsde
+c
+ dgde=rkfac*(half*epptdt(6)+ae*sdev0(6)-f4*sdevt(6))
+ dsde=f1*(one-dgde*(c2*sdevt(6)+
+ & f3*(epptdt(6)-f4*sdevt(6)+ae*sdev0(6))))
+ dmat(iddmat(6,6))=half*dsde
+c
return
end
c
@@ -633,10 +772,8 @@
c
c... local variables
c
- double precision ptmp(20)
+ double precision ptmp(10)
c
-cdebug write(6,*) "Hello from prestr_mat_7_f!"
-c
call dcopy(nprop,prop,ione,ptmp,ione)
if(ipauto.eq.ione) then
ptmp(2)=tyoungs
@@ -667,6 +804,7 @@
c
call dcopy(nstate,state,ione,sout,ione)
call dcopy(nstate,dstate,ione,sout(nstatesmax+ione),ione)
+c
return
end
c
@@ -697,9 +835,9 @@
double precision sub
data sub/-1.0d0/
c
-c... local variables
+cdebug write(6,*) "Hello from update_state_7_f!"
c
- call daxpy(nstate,sub,state,ione,dstate,ione)
+ call daxpy(nstate-6,sub,state,ione,dstate,ione)
call daxpy(nstate,one,dstate,ione,state,ione)
c
return
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 2007-03-19 16:50:26 UTC (rev 6287)
+++ short/3D/PyLith/branches/pylith-0.8/pylith3d/libpylith3d/mat_8.f 2007-03-19 16:56:29 UTC (rev 6288)
@@ -43,6 +43,8 @@
& ierr,errstrng)
c
c... subroutine to output material properties for material model 8.
+c This model is based on an ESF formulation although an actual
+c Effective Stress Function is not needed for a linear model.
c
c Error codes:
c 4: Write error
@@ -105,7 +107,7 @@
return
c
700 format(/,5x,"Material number: ",i7,/,5x,
- & "Material type: ",a37,/,5x,
+ & "Material type: ",a41,/,5x,
& "Number of properties: ",i7,/)
710 format(15x,a15,3x,1pe15.8)
c
Modified: short/3D/PyLith/branches/pylith-0.8/pylith3d/libpylith3d/mat_9.f
===================================================================
--- short/3D/PyLith/branches/pylith-0.8/pylith3d/libpylith3d/mat_9.f 2007-03-19 16:50:26 UTC (rev 6287)
+++ short/3D/PyLith/branches/pylith-0.8/pylith3d/libpylith3d/mat_9.f 2007-03-19 16:56:29 UTC (rev 6288)
@@ -30,7 +30,7 @@
c... Program segment to define material model 9:
c
c Model number: 9
-c Model name: IsotropicPowerLawMaxwellViscoelasticESF
+c Model name: IsotropicPowerLawMaxwellViscoelasticZT
c Number material properties: 5
c Number state variables: 18
c Tangent matrix varies with state: True
@@ -62,13 +62,14 @@
c
c... local constants
c
- character labelp(5)*21,modelname*40
+ character labelp(5)*21,modelname*44
data labelp/"Density",
& "Young's modulus",
& "Poisson's ratio",
& "Power-law exponent",
& "Viscosity coefficient"/
- parameter(modelname="Isotropic Powerlaw Maxwell Viscoelastic ESF")
+ parameter(modelname=
+ & "Isotropic Power-Law Maxwell Viscoelastic Z&T")
integer mattype
parameter(mattype=9)
c
@@ -107,7 +108,7 @@
return
c
700 format(/,5x,"Material number: ",i7,/,5x,
- & "Material type: ",a40,/,5x,
+ & "Material type: ",a44,/,5x,
& "Number of properties: ",i7,/)
710 format(15x,a21,3x,1pe15.8)
c
@@ -241,6 +242,7 @@
c... parameter definitions
c
include "ndimens.inc"
+ include "nconsts.inc"
include "rconsts.inc"
c
c... subroutine arguments
@@ -260,8 +262,9 @@
c
c... local variables
c
- double precision e,pr,anpwr,emhu,f1,f2,gam,ae,rmu
+ double precision e,pr,anpwr,emhu,rmu
double precision sdev(nstr),sinv1,steff
+ double precision cmat(nddmat)
c
c... included variable definitions
c
@@ -271,40 +274,63 @@
c
cdebug write(6,*) "Hello from td_matinit_9_f!"
c
- call fill(dmat,zero,nddmat)
- e=prop(2)
- pr=prop(3)
- anpwr=prop(4)
- emhu=prop(5)
- ae=(one+pr)/e
- rmu=half*e/(one+pr)
- f1=third*e/(one-two*pr)
- call invar(sdev,sinv1,steff,state)
- gam=half*((steff/emhu)**(anpwr-one))/emhu
+c
+c... get stress invariants and compute elastic material matrix
+c
+ if(iopt.eq.1) then
+ call invar(sdev,sinv1,steff,state)
+ else
+ call invar(sdev,sinv1,steff,dstate)
+ end if
+ call elas_mat_9(dmat,prop,iddmat,nprop,ierr,errstrng)
tmax=big
- if(steff.ne.zero) tmax=((emhu/steff)**(anpwr-one))*emhu/rmu
- f2=third/(ae+deltp*gam)
- 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))
+c
+c... if second deviatoric invariant is zero, use only elastic solution
+c
+ if(steff.eq.zero) then
+ return
+c
+c... otherwise, invert elastic matrix and augment it by the viscous
+c Jacobian.
+c
+ else
+c
+c... invert elastic matrix
+c
+ call invsymp(nddmat,nstr,dmat,ierr,errstrng)
+ if(ierr.ne.izero) return
+c
+c... define material properties
+c
+ e=prop(2)
+ pr=prop(3)
+ anpwr=prop(4)
+ emhu=prop(5)
+ rmu=half*e/(one+pr)
+c
+c... compute viscous Jacobian and add it to inverted elastic matrix
+c
+ call dbds_9(cmat,nddmat,iddmat,sdev,nstr,steff,anpwr,emhu,
+ & alfap,deltp)
+ call daxpy(nddmat,one,cmat,ione,dmat,ione)
+c
+c... invert augmented matrix
+c
+ call invsymp(nddmat,nstr,dmat,ierr,errstrng)
+ if(ierr.ne.izero) return
+ tmax=((emhu/steff)**(anpwr-one))*emhu/rmu
+ end if
return
end
c
c
- subroutine td_strs_9(state,dstate,state0,ee,scur,dmat,prop,
- & rtimdat,rgiter,iter,ntimdat,iddmat,tmax,nstate,nstate0,nprop,
- & matchg,ierr,errstrng)
+ subroutine jaccmp_9(scur,n,fvec,nfjsize,fjac,rpar,nrpar,ipar,
+ & nipar,ierr,errstrng)
c
-c... subroutine to compute the current values for stress, total strain,
-c and viscous strain increment 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... subroutine to form the Jacobian used in finding the zero of the
+c stress function.
+c** Note that in this version fjac is always assumed to be a symmetric
+c matrix stored in packed format.
c
include "implicit.inc"
c
@@ -313,252 +339,341 @@
include "ndimens.inc"
include "nconsts.inc"
include "rconsts.inc"
- integer nrpar,nipar
- parameter(nrpar=9,nipar=1)
+ include "parmat_9.inc"
c
c... subroutine arguments
c
- integer iter,nstate,nstate0,nprop,ierr
- integer iddmat(nstr,nstr)
+ integer n,nfjsize,nrpar,nipar,ierr
+c*** note that the dimension below is a bit kludgy and does not allow
+c*** anything else to be put in the ipar array
+ integer ipar(nstr,nstr)
+ double precision scur(n),fvec(n),fjac(nfjsize),rpar(nrpar)
character errstrng*(*)
- logical matchg
- double precision state(nstate),dstate(nstate),state0(nstate0)
- double precision ee(nstr),scur(nstr),dmat(nddmat),prop(nprop),tmax
c
-c... included dimension and type statements
+c... local variables
c
- include "rtimdat_dim.inc"
- include "rgiter_dim.inc"
- include "ntimdat_dim.inc"
+ integer i,j
+ctest integer iddmat(nstr,nstr)
+ctest equivalence(ipar(indiddmat),iddmat)
+ double precision cmat(nddmat)
+ double precision sdev(nstr),sinv1,steff
c
-c... external functions
+cdebug write(6,*) "Hello from jaccmp_9_f!"
c
- double precision sprod
- external sprod
c
-c... local constants
+c... get stress invariants and a copy of inverted elasticity matrix
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/
+ call invar(sdev,sinv1,steff,scur)
+ call dcopy(nddmat,rpar(inddmati),ione,fjac,ione)
c
+c... if second deviatoric invariant is zero, use only elastic solution
+c
+ if(steff.eq.zero) then
+ return
+c
+c... otherwise, augment inverted elastic matrix by the viscous Jacobian.
+c
+ else
+c
+c... compute viscous Jacobian and add it to inverted elastic matrix
+c
+ call dbds_9(cmat,nddmat,ipar,sdev,nstr,steff,
+ & rpar(indanpwr),rpar(indemhu),rpar(indalfap),rpar(inddeltp))
+ call daxpy(nddmat,one,cmat,ione,fjac,ione)
+ end if
+ return
+ end
+c
+c
+ subroutine dbds_9(cmat,nddmat,iddmat,sdev,nstr,steff,anpwr,emhu,
+ & alfap,deltp)
+c
+c... subroutine to compute derivatives of viscous strain rate with
+c respect to stress.
+c
+ implicit none
+c
+c... parameter definitions
+c
+ include "nconsts.inc"
+ include "rconsts.inc"
+c
+c... subroutine arguments
+c
+ integer nddmat,nstr
+ integer iddmat(nstr,nstr)
+ double precision cmat(nddmat),sdev(nstr),steff,anpwr,emhu
+ double precision alfap,deltp
+c
c... local variables
c
- integer i
- integer ipar(nipar)
- double precision rpar(nrpar)
- double precision e,pr,anpwr,emhu,rmu,ae,tf
- double precision emeantdt,smeant,smeantdt,smean0
- double precision sinv1t,sefft,sinv10,seff0,sefftdt
- double precision sefftau,gamtau,f1,f2,sdevtdt,sdevtau
- double precision epptdt(6),sdevt(6),sdev0(6)
-
+ double precision sigma,sigmap,sxx,syy,szz,sxy,syz,sxz
c
-c... included variable definitions
+cdebug write(6,*) "Hello from dbds_9!"
c
- include "rtimdat_def.inc"
- include "rgiter_def.inc"
- include "ntimdat_def.inc"
+ call fill(cmat,zero,nddmat)
+ sigma=(steff/emhu)**(anpwr-one)
+ sigmap=(anpwr-one)*sigma
+ sigma=half*alfap*deltp*sigma/emhu
+ sigmap=fourth*alfap*deltp*sigmap/emhu
+ sxx=sdev(1)/steff
+ syy=sdev(2)/steff
+ szz=sdev(3)/steff
+ sxy=two*sdev(4)/steff
+ syz=two*sdev(5)/steff
+ sxz=two*sdev(6)/steff
c
-cdebug write(6,*) "Hello from td_strs_9!"
+c... compute viscous Jacobian
c
- ierr=0
+ cmat(iddmat(1,1))=sigma*two*third+sigmap*sxx*sxx
+ cmat(iddmat(1,2))=-sigma*third +sigmap*sxx*syy
+ cmat(iddmat(1,3))=-sigma*third +sigmap*sxx*szz
+ cmat(iddmat(1,4))= sigmap*sxx*sxy
+ cmat(iddmat(1,5))= sigmap*sxx*syz
+ cmat(iddmat(1,6))= sigmap*sxx*sxz
+ cmat(iddmat(2,2))=sigma*two*third+sigmap*syy*syy
+ cmat(iddmat(2,3))=-sigma*third +sigmap*syy*szz
+ cmat(iddmat(2,4))= sigmap*syy*sxy
+ cmat(iddmat(2,5))= sigmap*syy*syz
+ cmat(iddmat(2,6))= sigmap*syy*sxz
+ cmat(iddmat(3,3))=sigma*two*third+sigmap*szz*szz
+ cmat(iddmat(3,4))= sigmap*szz*sxy
+ cmat(iddmat(3,5))= sigmap*szz*syz
+ cmat(iddmat(3,6))= sigmap*szz*sxz
+ cmat(iddmat(4,4))=sigma*two +sigmap*sxy*sxy
+ cmat(iddmat(4,5))= sigmap*sxy*syz
+ cmat(iddmat(4,6))= sigmap*sxy*sxz
+ cmat(iddmat(5,5))=sigma*two +sigmap*syz*syz
+ cmat(iddmat(5,6))= sigmap*syz*sxz
+ cmat(iddmat(6,6))=sigma*two +sigmap*sxz*sxz
+ return
+ end
c
-c... define material properties
c
- e=prop(2)
- pr=prop(3)
- anpwr=prop(4)
- emhu=prop(5)
- rmu=half*e/(one+pr)
- ae=(one+pr)/e
- tf=deltp*(one-alfap)
+ subroutine betacmp_9(strs,beta,sdev,seff,sinv1,emhu,anpwr,nstr)
c
-c... define stress and strain quantities
+c... subroutine to compute viscous strain rate vector beta
+c Routine also returns stress invariants
c
- emeantdt=third*(ee(1)+ee(2)+ee(3))
- smeantdt=e*emeantdt/(one-two*pr)
- call invar(sdevt,sinv1t,sefft,state)
- smeant=third*sinv1t
- call invar(sdev0,sinv10,seff0,state0)
- smean0=third*sinv10
- do i=1,nstr
- epptdt(i)=eng(i)*ee(i)-diag(i)*emeantdt-state(i+12)
- end do
+ implicit none
c
-c... define parameters needed by effective stress function
+c... parameter definitions
c
- rpar(1)=ae
- rpar(2)=anpwr
- rpar(3)=emhu
- rpar(4)=sprod(epptdt,epptdt)
- rpar(5)=two*ae*sprod(epptdt,sdev0)+ae*ae*seff0*seff0
- rpar(6)=two*tf*(sprod(epptdt,sdevt)+ae*sprod(sdevt,sdev0))
- rpar(7)=sefft
- rpar(8)=deltp
- rpar(9)=alfap
+ include "nconsts.inc"
+ include "rconsts.inc"
c
-c... call effective stress driver routine
+c... subroutine arguments
c
- call esfcomp_9(sefftdt,stol,rpar,nrpar,ipar,nipar,ierr,errstrng)
- if(ierr.ne.izero) return
+ integer nstr
+ double precision strs(nstr),beta(nstr),sdev(nstr),seff,sinv1
+ double precision emhu,anpwr
c
-c... compute quantities at intermediate time tau used to compute
-c values at end of time step.
+c... local variables
c
- sefftau=(one-alfap)*sefft+alfap*sefftdt
- gamtau=half*((sefftau/emhu)**(anpwr-one))/emhu
- f1=one/(ae+alfap*deltp*gamtau)
- f2=tf*gamtau
- tmax=big
- if(sefftdt.ne.zero) tmax=((emhu/sefftdt)**(anpwr-one))*emhu/rmu
+ double precision rm
c
-c... compute new stresses and store stress and strain values in dstate
+c... compute stress invariants and constants for loop
c
- do i=1,nstr
- sdevtdt=f1*(epptdt(i)-f2*sdevt(i)+ae*sdev0(i))
- dstate(i)=sdevtdt+diag(i)*(smeantdt+smean0)
- scur(i)=dstate(i)
- sdevtau=(one-alfap)*sdevt(i)+alfap*sdevtdt
- dstate(i+6)=ee(i)
- dstate(i+12)=deltp*gamtau*sdevtau
- end do
+ call fill(beta,zero,nstr)
+ call invar(sdev,sinv1,seff,strs)
+ rm=half*(seff/emhu)**(anpwr-one)/emhu
c
+c... compute strain rate components
+c
+ call daxpy(nstr,rm,sdev,ione,beta,ione)
return
end
c
c
- subroutine esfcomp_9(sefftdt,stol,rpar,nrpar,ipar,nipar,
- & ierr,errstrng)
+ subroutine td_strs_9(state,dstate,state0,ee,scur,dmat,prop,
+ & rtimdat,rgiter,iter,ntimdat,iddmat,tmax,nstate,nstate0,nprop,
+ & matchg,ierr,errstrng)
c
-c... driver routine to compute the effective stress at the current
-c time step.
+c... subroutine to compute the current values for stress, total strain,
+c and viscous strain increment 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
include "implicit.inc"
c
c... parameter definitions
c
+ include "ndimens.inc"
+ include "nconsts.inc"
include "rconsts.inc"
+ include "parmat_9.inc"
+ integer nrpar,nipar
+ parameter(nrpar=31,nipar=36)
c
c... subroutine arguments
c
- integer nrpar,nipar,ierr
- integer ipar(nipar)
- double precision sefftdt,stol,rpar(nrpar)
+ integer iter,nstate,nstate0,nprop,ierr
+ integer iddmat(nstr,nstr)
character errstrng*(*)
+ logical matchg
+ double precision state(nstate),dstate(nstate),state0(nstate0)
+ double precision ee(nstr),scur(nstr),dmat(nddmat),prop(nprop),tmax
c
c... included dimension and type statements
c
+ include "rtimdat_dim.inc"
+ include "rgiter_dim.inc"
+ include "ntimdat_dim.inc"
c
-c... intrinsic functions
+c... external functions
c
- intrinsic sqrt,max
+ double precision sprod,dasum
+ external sprod,funcv_9,jaccmp_9,dasum
c
-c... external routines
+c... local constants
c
- external esf_9
- double precision rtsafe
- external rtsafe
+ double precision sub
+ data sub/-1.0d0/
c
c... local variables
c
- double precision sefft,sefflo,seffhi,xacc
+ integer i
+ double precision e,pr,anpwr,emhu,rmu
+ double precision test0,rmt,rmtdt
+ double precision sdev(nstr),betat(nstr),betatdt(nstr),sinv1,seff
+ double precision rpar(nrpar)
+ logical check
c
c... included variable definitions
c
+ include "rtimdat_def.inc"
+ include "rgiter_def.inc"
+ include "ntimdat_def.inc"
c
-cdebug write(6,*) "Hello from esfcomp_9!"
+cdebug write(6,*) "Hello from td_strs_9!"
c
ierr=0
+ test0=dasum(nstate0,state0,ione)
c
-c... use stress from previous step as initial high value, defaulting to
-c stol if this is zero.
+c... define material properties
c
- sefft=rpar(7)
- seffhi=max(sefft,stol)
- sefflo=half*seffhi
- xacc=max(stol*half*(seffhi-sefflo),stol)
+ e=prop(2)
+ pr=prop(3)
+ anpwr=prop(4)
+ emhu=prop(5)
+ rmu=half*e/(one+pr)
+ tmax=big
+ rmt=deltp*(one-alfap)
+ rmtdt=deltp*alfap
c
-c... bracket the root
+c... for first iteration, use stresses from previous step as initial
+c guess, otherwise use current estimate.
c
- call zbrac(esf_9,sefflo,seffhi,rpar,nrpar,ipar,nipar,
- & ierr,errstrng)
- if(ierr.ne.0) return
+ if(iter.eq.ione) call dcopy(nstr,state,ione,dstate,ione)
c
-c... compute effective stress
+c... compute constant part of iterative solution equation.
c
- sefftdt=rtsafe(esf_9,sefflo,seffhi,xacc,rpar,nrpar,ipar,nipar,
- & ierr,errstrng)
+c... current total strain
+ call dcopy(nstr,ee,ione,rpar(indconst),ione)
+ call dscal(nstr,sub,rpar(indconst),ione)
+c... inverted elasticity matrix times initial stresses
+ call elas_mat_9(rpar(inddmati),prop,iddmat,nprop,ierr,errstrng)
+ call invsymp(nddmat,nstr,rpar(inddmati),ierr,errstrng)
+ if(ierr.ne.izero) return
+ if(test0.ne.zero) call dspmv("u",nstr,sub,rpar(inddmati),state0,
+ & ione,one,rpar(indconst),ione)
+c... viscous strain from previous time step
+ call daxpy(nstr,one,state(13),ione,rpar(indconst),ione)
+c... contribution to viscous strain from stresses at beginning of step
+ call betacmp_9(state,betat,sdev,seff,sinv1,emhu,anpwr,nstr)
+ call daxpy(nstr,rmt,betat,ione,rpar(indconst),ione)
+c
+c... finish defining parameters for stress computation then call Newton
+c routine to compute stresses.
+c
+ rpar(indalfap)=alfap
+ rpar(inddeltp)=deltp
+ rpar(indemhu)=emhu
+ rpar(indanpwr)=anpwr
+ call newt(dstate,nstr,rpar,nrpar,iddmat,nipar,funcv_9,jaccmp_9,
+ & check,ierr,errstrng)
+ if(ierr.ne.izero) return
+c
+c... update state variables for current stress estimates.
+c
+ call betacmp_9(dstate,betatdt,sdev,seff,sinv1,emhu,anpwr,nstr)
+ call dcopy(nstr,ee,ione,dstate(7),ione)
+ call dcopy(nstr,state(13),ione,dstate(13),ione)
+ call daxpy(nstr,rmt,betat,ione,dstate(13),ione)
+ call daxpy(nstr,rmtdt,betatdt,ione,dstate(13),ione)
+ call dcopy(nstr,dstate,ione,scur,ione)
+c
+c... compute current Maxwell time
+c
+ call invar(sdev,sinv1,seff,dstate)
+ if(seff.ne.zero) tmax=((emhu/seff)**(anpwr-one))*emhu/rmu
+c
return
end
c
c
- subroutine esf_9(seffcur,f,df,rpar,nrpar,ipar,nipar)
+ subroutine funcv_9(n,scur,r,rpar,nrpar,ipar,nipar,
+ & ierr,errstrng)
c
-c... subroutine to compute the effective stress function and it's
-c derivative for a given effective stress for a Maxwell power-law
-c viscoelastic material.
+c... routine to compute the vector of functions(r) evaluated at the
+c given stress state (contained in scur).
c
- include "implicit.inc"
+c... The rpar array contains:
+c 1-6: Constant part of vector, which includes current total
+c strain, viscous strain from previous step, stress-related
+c quantities from prefious time step, and prestress term.
+c 7-27: Inverted elastic material matrix.
+c 28: Integration parameter alfap.
+c 29: Time step size deltp.
+c 30: Viscosity coefficient emhu.
+c 31: Power-law exponent anpwr.
c
+c... The ipar array contains:
+c 1-36: The iddmat index array
+c
+ implicit none
+c
c... parameter definitions
c
+ include "ndimens.inc"
+ include "nconsts.inc"
include "rconsts.inc"
+ include "parmat_9.inc"
c
c... subroutine arguments
c
- integer nrpar,nipar
+ integer n,nrpar,nipar,ierr
integer ipar(nipar)
- double precision seffcur,f,df,rpar(nrpar)
+ double precision scur(n),r(n),rpar(nrpar)
+ character errstrng*(*)
c
-c... included dimension and type statements
+c... local data
c
+ double precision sub
+ data sub/-1.0d0/
c
c... local variables
c
- double precision ae,anpwr,emhu,a,b,c,sefft,deltp,alfap,d
- double precision f1,tf
- double precision gamtau,sefftau,dgam
+ double precision betatdt(nstr),rmtdt
+ double precision sdev(nstr),seff,sinv1
c
-c... included variable definitions
+c... viscous strain rate multiplier
c
+ rmtdt=rpar(indalfap)*rpar(inddeltp)
c
-cdebug write(6,*) "Hello from esf_9!"
+c... compute current viscous strain rate
c
+ call betacmp_9(scur,betatdt,sdev,seff,sinv1,
+ & rpar(indemhu),rpar(indanpwr),nstr)
c
-c... values from parameter list and a few other constants
+c... compute contributions to r and accumulate
c
- ae=rpar(1)
- anpwr=rpar(2)
- emhu=rpar(3)
- b=rpar(4)+rpar(5)
- c=rpar(6)
- sefft=rpar(7)
- deltp=rpar(8)
- alfap=rpar(9)
- f1=one-alfap
- tf=deltp*f1
- d=tf*sefft
+ call dcopy(nstr,rpar(indconst),ione,r,ione)
+ call daxpy(nstr,rmtdt,betatdt,ione,r,ione)
+ call dspmv("u",nstr,one,rpar(inddmati),scur,ione,one,r,ione)
c
-c... additional values dependent on current effective stress
-c
- sefftau=f1*sefft+alfap*seffcur
- gamtau=half*((sefftau/emhu)**(anpwr-one))/emhu
- a=ae+alfap*deltp*gamtau
-c
-c... effective stress function and it's derivative
-c
- f=a*a*seffcur*seffcur-b+c*gamtau-d*d*gamtau*gamtau
- dgam=half*alfap*(anpwr-one)*((sefftau/emhu)**(anpwr-two))/
- & (emhu*emhu)
- df=two*a*a*seffcur+(two*a*alfap*deltp*seffcur*seffcur+
- & c-two*d*d*gamtau)*dgam
-cdebug write(6,*) seffcur,f,df
-c
return
end
-
c
c
subroutine td_strs_mat_9(state,dstate,state0,ee,scur,dmat,prop,
@@ -577,8 +692,6 @@
include "ndimens.inc"
include "nconsts.inc"
include "rconsts.inc"
- integer nrpar,nipar
- parameter(nrpar=9,nipar=1)
c
c... subroutine arguments
c
@@ -597,27 +710,13 @@
c
c... external routines
c
- double precision sprod
- external sprod
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
- integer ipar(nipar)
- double precision rpar(nrpar)
- double precision e,pr,anpwr,emhu,rmu,ae,c1,c2,c3,a,c,d,rk1,rk2
- double precision emeantdt,smeant,smeantdt,smean0
- double precision sinv1t,sefft,sinv10,seff0,sefftdt
- double precision sefftau,gamtau,sdevtdt,sdevtau
- double precision f1,f2,f3,f4,rkdnm,rkfac,dgde,dsde
- double precision epptdt(6),sdevt(6),sdev0(6)
+ integer iopt
c
c... included variable definitions
c
@@ -625,123 +724,21 @@
include "rgiter_def.inc"
include "ntimdat_def.inc"
c
-c... rather than call td_strs_9, we duplicate the functionality here,
-c since otherwise there would be an enormous amount of duplicate
-c computations.
-c
ierr=0
+ iopt=2
c
-c... define material properties
+c... compute current stress estimates for time-dependent solution
c
- e=prop(2)
- pr=prop(3)
- anpwr=prop(4)
- emhu=prop(5)
- rmu=half*e/(one+pr)
- ae=(one+pr)/e
- c1=one-alfap
- c2=deltp*c1
- c3=third*e/(one-two*pr)
+ call td_strs_9(state,dstate,state0,ee,scur,dmat,prop,
+ & rtimdat,rgiter,iter,ntimdat,iddmat,tmax,nstate,nstate0,nprop,
+ & matchg,ierr,errstrng)
c
-c... define stress and strain quantities
+c... compute tangent stiffness corresponding to new stress estimates
c
- emeantdt=ee(1)+ee(2)+ee(3)
- smeantdt=c3*emeantdt
- emeantdt=third*emeantdt
- call invar(sdevt,sinv1t,sefft,state)
- smeant=third*sinv1t
- call invar(sdev0,sinv10,seff0,state0)
- smean0=third*sinv10
- do i=1,nstr
- epptdt(i)=eng(i)*ee(i)-diag(i)*emeantdt-state(i+12)
- end do
+ call td_matinit_9(state,dstate,state0,dmat,prop,rtimdat,
+ & rgiter,iopt,ntimdat,iddmat,tmax,nstate,nstate0,nprop,matchg,ierr,
+ & errstrng)
c
-c... define parameters needed by effective stress function
-c
- rpar(1)=ae
- rpar(2)=anpwr
- rpar(3)=emhu
- rpar(4)=sprod(epptdt,epptdt)
- rpar(5)=two*ae*sprod(epptdt,sdev0)+ae*ae*seff0*seff0
- rpar(6)=two*c2*(sprod(epptdt,sdevt)+ae*sprod(sdevt,sdev0))
- rpar(7)=sefft
- rpar(8)=deltp
- rpar(9)=alfap
-c
-c... call effective stress driver routine
-c
- call esfcomp_9(sefftdt,stol,rpar,nrpar,ipar,nipar,ierr,errstrng)
- if(ierr.ne.izero) return
-c
-c... compute quantities at intermediate time tau used to compute
-c values at end of time step.
-c
- sefftau=(one-alfap)*sefft+alfap*sefftdt
- gamtau=half*((sefftau/emhu)**(anpwr-one))/emhu
- a=ae+alfap*deltp*gamtau
- f1=one/a
- f2=c2*gamtau
-c
-c... compute new stresses and store stress and strain values in dstate
-c
- do i=1,nstr
- sdevtdt=f1*(epptdt(i)-f2*sdevt(i)+ae*sdev0(i))
- dstate(i)=sdevtdt+diag(i)*(smeantdt+smean0)
- scur(i)=dstate(i)
- sdevtau=(one-alfap)*sdevt(i)+alfap*sdevtdt
- dstate(i+6)=ee(i)
- dstate(i+12)=deltp*gamtau*sdevtau
- end do
- tmax=big
- if(sefftdt.ne.zero) tmax=((emhu/sefftdt)**(anpwr-one))*emhu/rmu
-c
-c... define some constants for tangent computation
-c
- rk1=half*(anpwr-one)*((sefftau/emhu)**(anpwr-two))/(emhu*emhu)
- rk2=sefftau-f1*sefft
- c=rpar(6)
- d=c2*rpar(7)
- f3=alfap*deltp*f1
- f4=gamtau*c2
- rkdnm=two*a*rk2*(alfap*deltp*rk1*rk2+a)/(alfap*alfap)+
- & rk1*(c-two*d*d*gamtau)
- rkfac=rk1/rkdnm
-c
-c... compute tangent matrix
-c
- dgde=rkfac*(half*epptdt(1)+ae*sdev0(1)-f4*sdevt(1))
- dsde=f1*(one-dgde*(c2*sdevt(1)+
- & f3*(epptdt(1)-f4*sdevt(1)+ae*sdev0(1))))
- dmat(iddmat(1,1))=c3+two*third*dsde
- dmat(iddmat(1,2))=c3-third*dsde
- dmat(iddmat(1,3))=dmat(iddmat(1,2))
-c
- dgde=rkfac*(half*epptdt(2)+ae*sdev0(2)-f4*sdevt(2))
- dsde=f1*(one-dgde*(c2*sdevt(2)+
- & f3*(epptdt(2)-f4*sdevt(2)+ae*sdev0(2))))
- dmat(iddmat(2,2))=c3+two*third*dsde
- dmat(iddmat(2,3))=c3-third*dsde
-c
- dgde=rkfac*(half*epptdt(3)+ae*sdev0(3)-f4*sdevt(3))
- dsde=f1*(one-dgde*(c2*sdevt(3)+
- & f3*(epptdt(3)-f4*sdevt(3)+ae*sdev0(3))))
- dmat(iddmat(3,3))=c3+two*third*dsde
-c
- dgde=rkfac*(half*epptdt(4)+ae*sdev0(4)-f4*sdevt(4))
- dsde=f1*(one-dgde*(c2*sdevt(4)+
- & f3*(epptdt(4)-f4*sdevt(4)+ae*sdev0(4))))
- dmat(iddmat(4,4))=half*dsde
-c
- dgde=rkfac*(half*epptdt(5)+ae*sdev0(5)-f4*sdevt(5))
- dsde=f1*(one-dgde*(c2*sdevt(5)+
- & f3*(epptdt(5)-f4*sdevt(5)+ae*sdev0(5))))
- dmat(iddmat(5,5))=half*dsde
-c
- dgde=rkfac*(half*epptdt(6)+ae*sdev0(6)-f4*sdevt(6))
- dsde=f1*(one-dgde*(c2*sdevt(6)+
- & f3*(epptdt(6)-f4*sdevt(6)+ae*sdev0(6))))
- dmat(iddmat(6,6))=half*dsde
-c
return
end
c
@@ -837,7 +834,7 @@
c
cdebug write(6,*) "Hello from update_state_9_f!"
c
- call daxpy(nstate-6,sub,state,ione,dstate,ione)
+ call daxpy(nstate,sub,state,ione,dstate,ione)
call daxpy(nstate,one,dstate,ione,state,ione)
c
return
Modified: short/3D/PyLith/branches/pylith-0.8/pylith3d/libpylith3d/matmod_def.f
===================================================================
--- short/3D/PyLith/branches/pylith-0.8/pylith3d/libpylith3d/matmod_def.f 2007-03-19 16:50:26 UTC (rev 6287)
+++ short/3D/PyLith/branches/pylith-0.8/pylith3d/libpylith3d/matmod_def.f 2007-03-19 16:56:29 UTC (rev 6288)
@@ -116,22 +116,22 @@
infmatmod(5,5) = izero
infmatmod(6,5) = 6
c
-c... Definition for power-law Maxwell viscoelastic material
+c... Definition for generalized linear Maxwell viscoelastic material
c
infmatmod(1,6) = ione
- infmatmod(2,6) = 18
- infmatmod(3,6) = ifive
+ infmatmod(2,6) = 30
+ infmatmod(3,6) = 9
infmatmod(4,6) = ione
- infmatmod(5,6) = ione
+ infmatmod(5,6) = izero
infmatmod(6,6) = 6
c
-c... Definition for generalized linear Maxwell viscoelastic material
+c... Definition for power-law Maxwell viscoelastic material (ESF)
c
infmatmod(1,7) = ione
- infmatmod(2,7) = 30
- infmatmod(3,7) = 9
+ infmatmod(2,7) = 18
+ infmatmod(3,7) = ifive
infmatmod(4,7) = ione
- infmatmod(5,7) = izero
+ infmatmod(5,7) = ione
infmatmod(6,7) = 6
c
c... Definition for linear Maxwell viscoelastic material (ESF version)
@@ -143,7 +143,7 @@
infmatmod(5,8) = izero
infmatmod(6,8) = 6
c
-c... Definition for power-law Maxwell viscoelastic material (ESF)
+c... Definition for power-law Maxwell viscoelastic material (Z&T)
c
infmatmod(1,9) = ione
infmatmod(2,9) = 18
Modified: short/3D/PyLith/branches/pylith-0.8/pylith3d/libpylith3d/newt.f
===================================================================
--- short/3D/PyLith/branches/pylith-0.8/pylith3d/libpylith3d/newt.f 2007-03-19 16:50:26 UTC (rev 6287)
+++ short/3D/PyLith/branches/pylith-0.8/pylith3d/libpylith3d/newt.f 2007-03-19 16:56:29 UTC (rev 6288)
@@ -36,7 +36,7 @@
c... local variables
c
double precision fvec(NP)
-CU USES fdjac,fmin,lnsrch,lubksb,ludcmp
+CU USES fdjac,fmin,lnsrch
INTEGER i,its,j,indx(NP)
double precision d,den,f,fold,stpmax,sum,temp,test,fjac(nddmat)
double precision g(NP),p(NP),xold(NP),fmin
More information about the cig-commits
mailing list