[cig-commits] r5779 -
short/3D/PyLith/branches/pylith-0.8/pylith3d/libpylith3d
willic3 at geodynamics.org
willic3 at geodynamics.org
Fri Jan 12 13:19:06 PST 2007
Author: willic3
Date: 2007-01-12 13:19:06 -0800 (Fri, 12 Jan 2007)
New Revision: 5779
Modified:
short/3D/PyLith/branches/pylith-0.8/pylith3d/libpylith3d/mat_6.f
Log:
New implementation of Maxwell power-law material model, using
formulation based on Zienkiewicz & Taylor rather than original
ESF-based version. Routines compile but have not yet been tested.
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-01-12 21:17:33 UTC (rev 5778)
+++ short/3D/PyLith/branches/pylith-0.8/pylith3d/libpylith3d/mat_6.f 2007-01-12 21:19:06 UTC (rev 5779)
@@ -32,7 +32,7 @@
c Model number: 6
c Model name: IsotropicPowerLawMaxwellViscoelastic
c Number material properties: 5
-c Number state variables: 3
+c Number state variables: 18
c Tangent matrix varies with state: True
c Material properties: Density
c Young's modulus
@@ -224,7 +224,7 @@
c
c
subroutine td_matinit_6(state,dstate,state0,dmat,prop,rtimdat,
- & rgiter,ntimdat,iddmat,tmax,nstate,nstate0,nprop,matchg,ierr,
+ & rgiter,iopt,ntimdat,iddmat,tmax,nstate,nstate0,nprop,matchg,ierr,
& errstrng)
c
c... subroutine to form the material matrix for an integration point
@@ -241,11 +241,12 @@
c... parameter definitions
c
include "ndimens.inc"
+ include "nconsts.inc"
include "rconsts.inc"
c
c... subroutine arguments
c
- integer nstate,nstate0,nprop,ierr
+ integer iopt,nstate,nstate0,nprop,ierr
integer iddmat(nstr,nstr)
character errstrng*(*)
double precision state(nstate),dstate(nstate)
@@ -260,8 +261,10 @@
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 sigma,sigmap,sxx,syy,szz,sxy,syz,sxz
+ double precision cmat(nddmat)
c
c... included variable definitions
c
@@ -271,35 +274,121 @@
c
cdebug write(6,*) "Hello from td_matinit_6_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_6(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
+ call dpptrf('u',nstr,dmat,ierr)
+ if(ierr.ne.izero) then
+ errstrng="td_matinit_6(1)"
+ if(ierr.lt.izero) then
+ ierr=117
+ else
+ ierr=118
+ end if
+ end if
+ call dpptri('u',nstr,dmat,ierr)
+ if(ierr.ne.izero) then
+ errstrng="td_matinit_6(2)"
+ if(ierr.lt.izero) then
+ ierr=117
+ else
+ ierr=119
+ end if
+ end if
+ call fill(cmat,zero,nddmat)
+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... get stress invariants and compute factors for Jacobian
+c
+ 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
+c
+c... add this to inverted elastic matrix and invert again
+c
+ call daxpy(nddmat,one,cmat,ione,dmat,ione)
+ call dpptrf('u',nstr,dmat,ierr)
+ if(ierr.ne.izero) then
+ errstrng="td_matinit_6(3)"
+ if(ierr.lt.izero) then
+ ierr=117
+ else
+ ierr=118
+ end if
+ end if
+ call dpptri('u',nstr,dmat,ierr)
+ if(ierr.ne.izero) then
+ errstrng="td_matinit_6(4)"
+ if(ierr.lt.izero) then
+ ierr=117
+ else
+ ierr=119
+ end if
+ end if
+ tmax=((emhu/steff)**(anpwr-one))*emhu/rmu
+ end if
return
end
c
c
subroutine td_strs_6(state,dstate,state0,ee,scur,dmat,prop,
- & rtimdat,rgiter,ntimdat,iddmat,tmax,nstate,nstate0,nprop,matchg,
- & ierr,errstrng)
+ & 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.
@@ -313,12 +402,10 @@
include "ndimens.inc"
include "nconsts.inc"
include "rconsts.inc"
- integer nrpar,nipar
- parameter(nrpar=9,nipar=1)
c
c... subroutine arguments
c
- integer nstate,nstate0,nprop,ierr
+ integer iter,nstate,nstate0,nprop,ierr
integer iddmat(nstr,nstr)
character errstrng*(*)
logical matchg
@@ -338,21 +425,15 @@
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,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 e,pr,anpwr,emhu,rmu
+ double precision sinv1t,sefft,sinv1tdt,sefftdt,beta,betat,betatdt
+ double precision rmt,rmtdt
+ double precision sdevt(nstr),sdevtdt(nstr),eetdt(nstr)
+ double precision dtmp(nddmat)
c
c... included variable definitions
@@ -372,198 +453,52 @@
anpwr=prop(4)
emhu=prop(5)
rmu=half*e/(one+pr)
- ae=(one+pr)/e
- tf=deltp*(one-alfap)
+ tmax=big
c
-c... define stress and strain quantities
+c... for first iteration, beta is computed using only the stresses at
+c time t. Otherwise, the stress estimates at time t + dt are also
+c used.
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
+ if(iter.eq.ione) then
+ do i=1,nstr
+ beta=half*(sefft/emhu)**(anpwr-one)*sdevt(i)/emhu
+ dstate(i+12)=state(i+12)+deltp*beta
+ eetdt(i)=ee(i)-dstate(i+12)
+ end do
+ else
+ rmt=one-alfap
+ rmtdt=alfap
+ call invar(sdevtdt,sinv1tdt,sefftdt,dstate)
+ do i=1,nstr
+ betat=half*(sefft/emhu)**(anpwr-one)*sdevt(i)/emhu
+ betatdt=half*(sefftdt/emhu)**(anpwr-one)*sdevtdt(i)/emhu
+ dstate(i+12)=state(i+12)+deltp*(rmt*betat+rmtdt*betatdt)
+ eetdt(i)=ee(i)-dstate(i+12)
+ end do
+ end if
c
-c... define parameters needed by effective stress function
+c... compute current stress estimates and copy state variables into the
+c appropriate spot
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
+ call elas_mat_6(dtmp,prop,iddmat,nprop,ierr,errstrng)
+ call dcopy(nstr,ee,ione,dstate(7),ione)
+ call dcopy(nstr,state0,ione,dstate,ione)
+ call dspmv("u",nstr,one,dtmp,eetdt,ione,one,dstate,ione)
+ call dcopy(nstr,dstate,ione,scur,ione)
c
-c... call effective stress driver routine
+c... compute current Maxwell time
c
- call esfcomp_6(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
+ call invar(sdevtdt,sinv1tdt,sefftdt,dstate)
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
- 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
-c
return
end
c
c
- subroutine esfcomp_6(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_6
- 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_6!"
-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_6,sefflo,seffhi,rpar,nrpar,ipar,nipar,
- & ierr,errstrng)
- if(ierr.ne.0) return
-c
-c... compute effective stress
-c
- sefftdt=rtsafe(esf_6,sefflo,seffhi,xacc,rpar,nrpar,ipar,nipar,
- & ierr,errstrng)
- return
- end
-c
-c
- subroutine esf_6(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_6!"
-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_6(state,dstate,state0,ee,scur,dmat,prop,
- & rtimdat,rgiter,ntimdat,iddmat,tmax,nstate,nstate0,nprop,matchg,
- & ierr,errstrng)
+ & rtimdat,rgiter,iter,ntimdat,iddmat,tmax,nstate,nstate0,nprop,
+ & matchg,ierr,errstrng)
c
c... subroutine to compute the current stress and updated material
c matrix for the time-dependent solution.
@@ -577,12 +512,10 @@
include "ndimens.inc"
include "nconsts.inc"
include "rconsts.inc"
- integer nrpar,nipar
- parameter(nrpar=9,nipar=1)
c
c... subroutine arguments
c
- integer nstate,nstate0,nprop,ierr
+ integer iter,nstate,nstate0,nprop,ierr
integer iddmat(nstr,nstr)
character errstrng*(*)
logical matchg
@@ -597,27 +530,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 +544,21 @@
include "rgiter_def.inc"
include "ntimdat_def.inc"
c
-c... rather than call td_strs_6, 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_6(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_6(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_6(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
More information about the cig-commits
mailing list