[cig-commits] r4778 -
short/3D/PyLith/branches/pylith-0.8/pylith3d/libpylith3d
willic3 at geodynamics.org
willic3 at geodynamics.org
Tue Oct 10 13:53:25 PDT 2006
Author: willic3
Date: 2006-10-10 13:53:24 -0700 (Tue, 10 Oct 2006)
New Revision: 4778
Modified:
short/3D/PyLith/branches/pylith-0.8/pylith3d/libpylith3d/mat_6.f
Log:
Started updating these routines for power-law Maxwell viscoelastic
rheology. Additional changes will be needed. Also, all other
material models will need to be updated to the new calling
conventions.
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 2006-10-10 20:37:07 UTC (rev 4777)
+++ short/3D/PyLith/branches/pylith-0.8/pylith3d/libpylith3d/mat_6.f 2006-10-10 20:53:24 UTC (rev 4778)
@@ -30,13 +30,15 @@
c... Program segment to define material model 6:
c
c Model number: 6
-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: IsotropicPowerLawMaxwellViscoelastic
+c Number material properties: 5
+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 Power-law exponent
+c Viscosity coefficient
c
subroutine mat_prt_6(prop,nprop,matnum,idout,idsk,kw,kp,
& ierr,errstrng)
@@ -45,7 +47,6 @@
c
c Error codes:
c 4: Write error
-c 101: Attempt to use undefined material model
c
include "implicit.inc"
c
@@ -61,20 +62,55 @@
c
c... local constants
c
- character labelp(3)*15,modelname*24
+ character labelp(5)*21,modelname*40
data labelp/"Density",
& "Young's modulus",
- & "Poisson's ratio"/
- parameter(modelname="???????")
+ & "Poisson's ratio",
+ & "Power-law exponent",
+ & "Viscosity coefficient"/
+ parameter(modelname="Isotropic Power-Law Maxwell Viscoelastic")
integer mattype
parameter(mattype=6)
c
-c... return error code, as this material is not yet defined
+c... local variables
c
- ierr=101
- errstrng="mat_prt_6"
+ integer i
c
+cdebug write(6,*) "Hello from mat_prt_6_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_6"
+ return
+c
+ 700 format(/,5x,"Material number: ",i7,/,5x,
+ & "Material type: ",a40,/,5x,
+ & "Number of properties: ",i7,/)
+ 710 format(15x,a21,3x,1pe15.8)
+c
end
c
c
@@ -100,21 +136,51 @@
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_6"
+ integer i,j
+ double precision e,pr,pr1,pr2,pr3,fac,dd,od,ss
+c
+cdebug write(6,*) "Hello from elas_mat_6_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
c
- subroutine elas_strs_6(state,state0,ee,scur,dmat,nstate,nstate0,
- & ierr,errstrng)
+ 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.
+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
@@ -125,17 +191,37 @@
c
c... subroutine arguments
c
- integer nstate,nstate0,ierr
- double precision state(nstate),state0(nstate0),ee(nstr),scur(nstr)
+ integer nprop,nstate,nstate0,ierr
+ double precision prop(nprop),state(nstate),state0(nstate0)
+ double precision ee(nstr),scur(nstr)
double precision dmat(nddmat)
character errstrng*(*)
c
-c... return error code, as this material is not yet defined
+c... local variables
c
- ierr=101
- errstrng="elas_strs_6"
+ double precision e,pr,anpwr,emhu,rmu
+ double precision sdev(nstr),sinv1,steff
+c
+cdebug write(6,*) "Hello from elas_strs_6_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)
+ 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
return
end
+c**************** finish fixing from here *****************
c
c
subroutine td_matinit_6(state,dstate,state0,dmat,prop,rtimdat,
More information about the cig-commits
mailing list