[cig-commits] commit: add cylindrical and spherical ductile zones

[Mercurial]

changeset:   206:ff574e3a9184
user:        Sylvain Barbot <sbarbot at ntu.edu.sg>
date:        Thu Jul 25 10:58:19 2013 +0200
files:       src/viscoelastic3d.f90
description:
add cylindrical and spherical ductile zones


diff -r f7b5f1f16a60 -r ff574e3a9184 src/viscoelastic3d.f90
--- a/src/viscoelastic3d.f90	Thu Jun 06 13:08:04 2013 +0800
+++ b/src/viscoelastic3d.f90	Thu Jul 25 10:58:19 2013 +0200
@@ -227,7 +227,7 @@
        TYPE(WEAK_STRUCT), INTENT(IN), DIMENSION(nz) :: zones
        REAL*8, INTENT(IN) :: x1,x2,x3,beta
 
-       REAL*8 :: dg,x,y,z,L,W,D,strike,dip,LM
+       REAL*8 :: dg,x,y,z,L,W,D,strike,dip,LM,r,rho
        REAL*8 :: cstrike,sstrike,cdip,sdip, &
                  xr,yr,zr,x2r,Wp,Lp,Dp,x1s,x2s,x3s
        INTEGER :: i
@@ -276,9 +276,27 @@
           x2s= sstrike*x1 +cstrike*x2
           x3s= sdip   *x2r+cdip   *x3
 
-          dgammadot0=dgammadot0+omega((x1s-xr)/D,beta) &
-                               *omega((x2s-yr)/W,beta) &
-                               *omega((x3s-zr)/L,beta)*dg
+          IF (D.GE.0 .AND. W.GE.0 .AND. L.GE.0) THEN
+             dgammadot0=dgammadot0+omega((x1s-xr)/D,beta) &
+                                  *omega((x2s-yr)/W,beta) &
+                                  *omega((x3s-zr)/L,beta)*dg
+          ELSE IF (D.LT.0 .AND. W.LT.0 .AND. L.LT.0) THEN
+             ! spherical ductile zone
+             r=SQRT(((x1s-xr)/D)**2+((x2s-yr)/W)**2+((x3s-zr)/L)**2)
+             dgammadot0=dgammadot0+omega(r,beta)*dg
+          ELSE IF (D.LT.0 .AND. W.LT.0) THEN
+             ! cylindrical ductile zone
+             rho=SQRT(((x1s-xr)/D)**2+((x2s-yr)/W)**2)
+             dgammadot0=dgammadot0+omega(rho,beta)*omega((x3s-zr)/L,beta)*dg
+          ELSE IF (D.LT.0 .AND. L.LT.0) THEN
+             ! cylindrical ductile zone
+             rho=SQRT(((x1s-xr)/D)**2+((x3s-zr)/L)**2)
+             dgammadot0=dgammadot0+omega(rho,beta)*omega((x2s-yr)/W,beta)*dg
+          ELSE IF (W.LT.0 .AND. L.LT.0) THEN
+             ! cylindrical ductile zone
+             rho=SQRT(((x2s-yr)/W)**2+((x3s-zr)/L)**2)
+             dgammadot0=dgammadot0+omega(rho,beta)*omega((x1s-xr)/D,beta)*dg
+          ENDIF
        END DO
 
     END FUNCTION dgammadot0