[cig-commits] r5304 - in
short/3D/PyLith/branches/pylith-0.8/pylith3d/libpylith3d: . includes
willic3 at geodynamics.org
willic3 at geodynamics.org
Thu Nov 16 13:08:23 PST 2006
Author: willic3
Date: 2006-11-16 13:08:22 -0800 (Thu, 16 Nov 2006)
New Revision: 5304
Modified:
short/3D/PyLith/branches/pylith-0.8/pylith3d/libpylith3d/includes/materials.inc
short/3D/PyLith/branches/pylith-0.8/pylith3d/libpylith3d/iterate.F
short/3D/PyLith/branches/pylith-0.8/pylith3d/libpylith3d/mat_7.f
short/3D/PyLith/branches/pylith-0.8/pylith3d/libpylith3d/matmod_def.f
short/3D/PyLith/branches/pylith-0.8/pylith3d/libpylith3d/read_stateout.f
Log:
More work on generalized Maxwell model.
Still need to finish mat_7.f.
Modified: short/3D/PyLith/branches/pylith-0.8/pylith3d/libpylith3d/includes/materials.inc
===================================================================
--- short/3D/PyLith/branches/pylith-0.8/pylith3d/libpylith3d/includes/materials.inc 2006-11-16 20:36:02 UTC (rev 5303)
+++ short/3D/PyLith/branches/pylith-0.8/pylith3d/libpylith3d/includes/materials.inc 2006-11-16 21:08:22 UTC (rev 5304)
@@ -30,15 +30,14 @@
c... materials.inc: A header that defines the maximum number
c of material models and state variables.
c At present, the maximum number of material models is 20 and the
-c maximum number of state variables is 24 (6 components each of
-c stress, strain, viscous strain, and plastic strain).
+c maximum number of state variables is 30.
c
c This header should be included before any subroutine arguments,
c as the parameters will generally be used to dimension the
c arguments.
c
integer nmatmodmax,nstatesmax
- parameter(nmatmodmax=20,nstatesmax=24)
+ parameter(nmatmodmax=20,nstatesmax=30)
c
c version
c $Id: materials.inc,v 1.2 2005/03/25 01:30:05 willic3 Exp $
Modified: short/3D/PyLith/branches/pylith-0.8/pylith3d/libpylith3d/iterate.F
===================================================================
--- short/3D/PyLith/branches/pylith-0.8/pylith3d/libpylith3d/iterate.F 2006-11-16 20:36:02 UTC (rev 5303)
+++ short/3D/PyLith/branches/pylith-0.8/pylith3d/libpylith3d/iterate.F 2006-11-16 21:08:22 UTC (rev 5304)
@@ -329,6 +329,7 @@
& skew,numrot, ! skew
& getshape,bmatrix, ! bbar
& ierr,errstrng) ! errcode
+ write(6,*) "tminmax: ",tminmax
else
call stress_drv(
& bintern,neq, ! force
@@ -343,6 +344,7 @@
& skew,numrot, ! skew
& getshape,bmatrix, ! bbar
& ierr,errstrng) ! errcode
+ write(6,*) "tminmax: ",tminmax
end if
if(ierr.ne.izero) return
c
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 2006-11-16 20:36:02 UTC (rev 5303)
+++ short/3D/PyLith/branches/pylith-0.8/pylith3d/libpylith3d/mat_7.f 2006-11-16 21:08:22 UTC (rev 5304)
@@ -32,7 +32,7 @@
c Model number: 7
c Model name: IsotropicLinearGenMaxwellViscoelastic
c Number material properties: 9
-c Number state variables: 15
+c Number state variables: 30
c Tangent matrix varies with state: True
c Material properties: Density
c Young's Modulus
@@ -91,14 +91,14 @@
integer i
double precision sfrac
c
-cdebug write(6,*) "Hello from mat_prt_5_f!"
+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_5"
+ errstrng="mat_prt_7"
return
end if
c
@@ -202,7 +202,12 @@
c
c state(1:6) = Cauchy stress
c state(7:12) = linear strain
-c state(13:15) = viscous state variable
+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
@@ -233,7 +238,7 @@
sfrac=prop(4)+prop(6)+prop(8)
if(sfrac.lt.0.0d0.or.sfrac.gt.1.0d0) then
ierr=116
- errstrng="mat_prt_5"
+ errstrng="elas_strs_7"
return
end if
c
@@ -364,15 +369,20 @@
sfrac=prop(4)+prop(6)+prop(8)
if(sfrac.lt.0.0d0.or.sfrac.gt.1.0d0) then
ierr=116
- errstrng="mat_prt_5"
+ errstrng="td_matinit_7"
return
end if
c
+c... get material properties
+c
call fill(dmat,zero,nddmat)
e=prop(2)
pr=prop(3)
rmu=half*e/(one+pr)
bulkm=third*e/(one-two*pr)
+c
+c... compute Prony series term
+c
frac0=one-sfrac
shfac=frac0
tmax=big
@@ -384,6 +394,9 @@
shfac=shfac+rmue*compdq(deltp,vise,rmue)
end if
end do
+c
+c... compute stress-strain matrix
+c
shfac=third*rmu*shfac
dmat(iddmat(1,1))=bulkm+two*shfac
dmat(iddmat(2,2))=dmat(iddmat(1,1))
@@ -430,17 +443,80 @@
include "rgiter_dim.inc"
include "ntimdat_dim.inc"
c
+c... external functions
+c
+ double precision compdq
+ external compdq
+c
+c... local variables
+c
+ integer i
+ double precision sfrac,e,pr,rmu,bulkm,rmue,vise,frac0,shfac
+ double precision etracetdt,emeantddt,smeantdt,sdevtdt,dq(3),qt(3)
+ double precision rtime(3)
+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_7_f!"
c
- ierr=101
- errstrng="td_strs_7"
+ ierr=0
+ sfrac=prop(4)+prop(6)+prop(8)
+ if(sfrac.lt.0.0d0.or.sfrac.gt.1.0d0) then
+ ierr=116
+ errstrng="td_strs_5"
+ return
+ end if
c
+c... define material properties and dilatational stress
+c
+ e=prop(2)
+ pr=prop(3)
+ rmu=half*e/(one+pr)
+ bulkm=third*e/(one-two*pr)
+ etracetdt=ee(1)+ee(2)+ee(3)
+ emeantdt=third*etracetdt
+ smeantdt=bulkm*(ee(1)+ee(2)+ee(3))
+ emeant=third*(state(7)+state(8)+state(9))
+c
+c... compute Prony series terms
+c
+ frac0=one-sfrac
+ shfac=frac0
+ tmax=big
+ do i=1,3
+ rmue=prop(2*(i-1)+4)
+ vise=prop(2*(i-1)+5)
+ rtime(i)=zero
+ if(rmue.ne.zero) then
+ rtime(i)=vise/rmue
+ tmax=min(tmax,rtime(i))
+ dstate(i+12)=qt(i)+compdq(deltp,vise,rmue)
+ qtdt(i)=rmue*dstate(i+12)
+ end if
+ end do
+ shfac=two*rmu*shfac
+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)
+ vise=prop(2*(j-1)+5)
+ if(rmue.ne.zero) then
+ qtdt=state(12+nstr*(j-1)+i)*exp(-deltp/rtime(j))+
+ & dq(j)*deltae
+ qt(i)=state(i+12)*exp(-deltp/rtime)
+ dstate(i+12)=qt(i)+compdq(deltp,vise,rmue)
+ qtdt(i)=rmue*dstate(i+12)
+c
return
end
c
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 2006-11-16 20:36:02 UTC (rev 5303)
+++ short/3D/PyLith/branches/pylith-0.8/pylith3d/libpylith3d/matmod_def.f 2006-11-16 21:08:22 UTC (rev 5304)
@@ -128,7 +128,7 @@
c... Definition for generalized linear Maxwell viscoelastic material
c
infmatmod(1,7) = ione
- infmatmod(2,7) = 15
+ infmatmod(2,7) = 30
infmatmod(3,7) = 9
infmatmod(4,7) = ione
infmatmod(5,7) = izero
Modified: short/3D/PyLith/branches/pylith-0.8/pylith3d/libpylith3d/read_stateout.f
===================================================================
--- short/3D/PyLith/branches/pylith-0.8/pylith3d/libpylith3d/read_stateout.f 2006-11-16 20:36:02 UTC (rev 5303)
+++ short/3D/PyLith/branches/pylith-0.8/pylith3d/libpylith3d/read_stateout.f 2006-11-16 21:08:22 UTC (rev 5304)
@@ -36,7 +36,7 @@
c desired for both the elastic and time-dependent solutions.
c
c The istatout array specifies output options for each individual
-c state variable. At present there are a maximum of 24 possible
+c state variable. At present there are a maximum of 30 possible
c state variables, and this number may increase with the addition
c of new material models. There are three types of state variable
c output:
More information about the cig-commits
mailing list