[cig-commits] commit: export fault stress to .xy format for gmt, add strike and dip components of stress for vtk output of fault stress.
Mercurial
hg at geodynamics.org
Sat Nov 26 15:00:16 PST 2011
changeset: 49:cf00c381e48f
tag: tip
user: Sylvain Barbot <sbarbot at caltech.edu>
date: Sat Nov 26 15:00:08 2011 -0800
files: export.f90 relax.f90
description:
export fault stress to .xy format for gmt, add strike and dip components of stress for vtk output of fault stress.
diff -r bf1575daaa90 -r cf00c381e48f export.f90
--- a/export.f90 Sat Nov 26 13:32:16 2011 -0800
+++ b/export.f90 Sat Nov 26 15:00:08 2011 -0800
@@ -1270,7 +1270,7 @@ END SUBROUTINE exportcreep
!!
!! \author sylvain barbot 06/06/11 - original form
!------------------------------------------------------------------
- SUBROUTINE exportvtk_rfaults_stress_init(sig,sx1,sx2,sx3,dx1,dx2,dx3, &
+ SUBROUTINE export_rfaults_stress_init(sig,sx1,sx2,sx3,dx1,dx2,dx3, &
nsop,sop)
INTEGER, INTENT(IN) :: sx1,sx2,sx3,nsop
TYPE(TENSOR), INTENT(IN), DIMENSION(sx1,sx2,sx3) :: sig
@@ -1300,7 +1300,144 @@ END SUBROUTINE exportcreep
END DO
- END SUBROUTINE exportvtk_rfaults_stress_init
+ END SUBROUTINE export_rfaults_stress_init
+
+ !------------------------------------------------------------------
+ !> subroutine ExportGMT_RFaults_Stress
+ !! creates a .vtp file (in the VTK PolyData XML format) containing
+ !! the rectangular faults. The faults are characterized with a set
+ !! of subsegments (rectangles) each associated with stress values.
+ !!
+ !! \author sylvain barbot 06/06/11 - original form
+ !------------------------------------------------------------------
+ SUBROUTINE exportgmt_rfaults_stress(sx1,sx2,sx3,dx1,dx2,dx3, &
+ nsop,sop,rffilename,convention,sig)
+ USE elastic3d
+ INTEGER, INTENT(IN) :: sx1,sx2,sx3,nsop
+ REAL*8, INTENT(IN) :: dx1,dx2,dx3
+ TYPE(SEGMENT_STRUCT), INTENT(INOUT), DIMENSION(nsop) :: sop
+ CHARACTER(80), INTENT(IN) :: rffilename
+ INTEGER, INTENT(IN), OPTIONAL :: convention
+ TYPE(TENSOR), INTENT(IN), DIMENSION(sx1,sx2,sx3), OPTIONAL :: sig
+
+ INTEGER :: iostatus,k,i1,i2,i3,conv
+ CHARACTER :: q
+ REAL*8 :: strike,dip,x1,x2,x3,cstrike,sstrike,cdip,sdip,L,W
+ ! segment normal vector, strike direction, dip direction
+ REAL*8, DIMENSION(3) :: n,s,d
+ ! local value of stress
+ TYPE(TENSOR) :: lsig
+ ! stress components
+ REAL*8 :: taun,taus,taustrike,taudip,taucoulomb
+ ! friction coefficient
+ REAL*8 :: friction
+ ! traction components
+ REAL*8, DIMENSION(3) :: t,ts
+
+ IF (0.GE.nsop) RETURN
+
+ ! double-quote character
+ q=char(34)
+
+ IF (PRESENT(convention)) THEN
+ conv=convention
+ ELSE
+ conv=0
+ END IF
+
+ OPEN (UNIT=15,FILE=rffilename,IOSTAT=iostatus,FORM="FORMATTED")
+ IF (iostatus>0) THEN
+ WRITE_DEBUG_INFO
+ PRINT '(a)', rffilename
+ STOP "could not open file for export"
+ END IF
+
+ DO k=1,nsop
+ ! friction coefficient
+ friction=sop(k)%friction
+
+ ! fault orientation
+ strike=sop(k)%strike
+ dip=sop(k)%dip
+
+ ! fault center position
+ x1=sop(k)%x
+ x2=sop(k)%y
+ x3=sop(k)%z
+
+ IF (PRESENT(sig)) THEN
+
+ CALL shiftedindex(x1,x2,x3,sx1,sx2,sx3,dx1,dx2,dx3,i1,i2,i3)
+ lsig=sig(i1,i2,i3)
+
+ IF (1.EQ.conv) THEN
+ lsig%s11=lsig%s11-sop(k)%sig0%s11
+ lsig%s12=lsig%s12-sop(k)%sig0%s12
+ lsig%s13=lsig%s13-sop(k)%sig0%s13
+ lsig%s22=lsig%s22-sop(k)%sig0%s22
+ lsig%s23=lsig%s23-sop(k)%sig0%s23
+ lsig%s33=lsig%s33-sop(k)%sig0%s33
+ END IF
+ ELSE
+ lsig=TENSOR(0.d0,0.d0,0.d0,0.d0,0.d0,0.d0)
+ END IF
+
+ ! fault dimension
+ W=sop(k)%width
+ L=sop(k)%length
+
+ cstrike=cos(strike)
+ sstrike=sin(strike)
+ cdip=cos(dip)
+ sdip=sin(dip)
+
+ ! surface normal vector components
+ n(1)=+cdip*cstrike
+ n(2)=-cdip*sstrike
+ n(3)=-sdip
+
+ ! strike-slip unit direction
+ s(1)=sstrike
+ s(2)=cstrike
+ s(3)=0._8
+
+ ! dip-slip unit direction
+ d(1)=+cstrike*sdip
+ d(2)=-sstrike*sdip
+ d(3)=+cdip
+
+ ! traction vector
+ t=lsig .tdot. n
+
+ ! signed normal component
+ taun=SUM(t*n)
+
+ ! shear traction
+ ts=t-taun*n
+
+ ! absolute value of shear component
+ taus=SQRT(SUM(ts*ts))
+
+ ! strike-direction shear component
+ taustrike=SUM(ts*s)
+
+ ! dip-direction shear component
+ taudip=SUM(ts*d)
+
+ ! Coulomb stress
+ taucoulomb=taus+friction*taun
+
+ WRITE (15,'("> -Z",5ES11.2)') taus, taun, taucoulomb, taustrike, taudip
+ WRITE (15,'(3ES11.2)') x1-d(1)*W/2-s(1)*L/2, x2-d(2)*W/2-s(2)*L/2
+ WRITE (15,'(3ES11.2)') x1-d(1)*W/2+s(1)*L/2, x2-d(2)*W/2+s(2)*L/2
+ WRITE (15,'(3ES11.2)') x1+d(1)*W/2+s(1)*L/2, x2+d(2)*W/2+s(2)*L/2
+ WRITE (15,'(3ES11.2)') x1+d(1)*W/2-s(1)*L/2, x2+d(2)*W/2-s(2)*L/2
+
+ END DO
+
+ CLOSE(15)
+
+ END SUBROUTINE exportgmt_rfaults_stress
!------------------------------------------------------------------
!> subroutine ExportVTK_RFaults_Stress
@@ -1439,10 +1576,10 @@ END SUBROUTINE exportcreep
" format=",a,"ascii",a,">")'),q,q,q,q,q,q,q,q
! fault edge coordinates
WRITE (15,'(12ES11.2)') &
- x1-d(1)*W/2-s(1)*L/2,x2-d(2)*W/2-s(2)*L/2,x3-d(3)*W/2-s(3)*L/2, &
- x1-d(1)*W/2+s(1)*L/2,x2-d(2)*W/2+s(2)*L/2,x3-d(3)*W/2+s(3)*L/2, &
- x1+d(1)*W/2+s(1)*L/2,x2+d(2)*W/2+s(2)*L/2,x3+d(3)*W/2+s(3)*L/2, &
- x1+d(1)*W/2-s(1)*L/2,x2+d(2)*W/2-s(2)*L/2,x3+d(3)*W/2-s(3)*L/2
+ x1-d(1)*W/2-s(1)*L/2, x2-d(2)*W/2-s(2)*L/2, x3-d(3)*W/2-s(3)*L/2, &
+ x1-d(1)*W/2+s(1)*L/2, x2-d(2)*W/2+s(2)*L/2, x3-d(3)*W/2+s(3)*L/2, &
+ x1+d(1)*W/2+s(1)*L/2, x2+d(2)*W/2+s(2)*L/2, x3+d(3)*W/2+s(3)*L/2, &
+ x1+d(1)*W/2-s(1)*L/2, x2+d(2)*W/2-s(2)*L/2, x3+d(3)*W/2-s(3)*L/2
WRITE (15,'(" </DataArray>")')
WRITE (15,'(" </Points>")')
@@ -1491,6 +1628,20 @@ END SUBROUTINE exportcreep
" NumberOfComponents=",a,"1",a, &
" format=",a,"ascii",a,">")'), q,q,q,q,q,q,q,q
WRITE (15,'(ES11.2)'), taucoulomb
+ WRITE (15,'(" </DataArray>")')
+
+ WRITE (15,'(" <DataArray type=",a,"Float32",a, &
+ " Name=",a,"stress in strike direction",a, &
+ " NumberOfComponents=",a,"1",a, &
+ " format=",a,"ascii",a,">")'), q,q,q,q,q,q,q,q
+ WRITE (15,'(ES11.2)'), taustrike
+ WRITE (15,'(" </DataArray>")')
+
+ WRITE (15,'(" <DataArray type=",a,"Float32",a, &
+ " Name=",a,"stress in dip direction",a, &
+ " NumberOfComponents=",a,"1",a, &
+ " format=",a,"ascii",a,">")'), q,q,q,q,q,q,q,q
+ WRITE (15,'(ES11.2)'), taudip
WRITE (15,'(" </DataArray>")')
WRITE (15,'(" </CellData>")')
diff -r bf1575daaa90 -r cf00c381e48f relax.f90
--- a/relax.f90 Sat Nov 26 13:32:16 2011 -0800
+++ b/relax.f90 Sat Nov 26 15:00:08 2011 -0800
@@ -211,8 +211,8 @@ PROGRAM relax
#ifdef VTK
CHARACTER(80) :: filename,title,name
CHARACTER(3) :: digit
+#endif
CHARACTER(4) :: digit4
-#endif
REAL*8 :: t,Dt,tm
! arrays
@@ -402,8 +402,11 @@ PROGRAM relax
in%x0,in%y0,in%lon0,in%lat0,in%zone,in%umult,in%wdir,0)
END IF
#endif
+ ! initialize stress conditions
+ CALL export_rfaults_stress_init(sig,in%sx1,in%sx2,in%sx3, &
+ in%dx1,in%dx2,in%dx3,in%nsop,in%sop)
+ WRITE (digit4,'(I4.4)') 0
#ifdef VTK
- WRITE (digit4,'(I4.4)') 0
IF (in%isoutputvtk .AND. in%isoutputstress) THEN
filename=trim(in%wdir)//"/sigma-"//digit4//".vtk"//char(0)
title="stress tensor field"//char(0)
@@ -411,16 +414,24 @@ PROGRAM relax
CALL exportvtk_tensors_legacy(sig,in%sx1,in%sx2,in%sx3/2,in%dx1,in%dx2,in%dx3, &
4,4,8,filename,title,name)
END IF
+ ! coseismic stress change on predefined planes for 3-D visualization w/ Paraview
filename=trim(in%wdir)//"/rfaults-sigma-"//digit4//".vtp"
CALL exportvtk_rfaults_stress(in%sx1,in%sx2,in%sx3,in%dx1,in%dx2,in%dx3, &
in%nsop,in%sop,filename,sig=sig)
- ! initialize stress conditions
- CALL exportvtk_rfaults_stress_init(sig,in%sx1,in%sx2,in%sx3, &
- in%dx1,in%dx2,in%dx3,in%nsop,in%sop)
+ ! postseismic stress change on predefined planes (zero by definition)
filename=trim(in%wdir)//"/rfaults-dsigma-"//digit4//".vtp"
CALL exportvtk_rfaults_stress(in%sx1,in%sx2,in%sx3,in%dx1,in%dx2,in%dx3, &
in%nsop,in%sop,filename)
#endif
+ ! coseismic stress change on predefined planes for gmt
+ filename=trim(in%wdir)//"/rfaults-sigma-"//digit4//".xy"
+ CALL exportgmt_rfaults_stress(in%sx1,in%sx2,in%sx3,in%dx1,in%dx2,in%dx3, &
+ in%nsop,in%sop,filename,sig=sig)
+ ! postseismic stress change on predefined planes for gmt (zero by definition)
+ filename=trim(in%wdir)//"/rfaults-dsigma-"//digit4//".xy"
+ CALL exportgmt_rfaults_stress(in%sx1,in%sx2,in%sx3,in%dx1,in%dx2,in%dx3, &
+ in%nsop,in%sop,filename)
+ ! time series of stress in ASCII format
CALL exportcoulombstress(sig,in%sx1,in%sx2,in%sx3,in%dx1,in%dx2,in%dx3, &
in%nsop,in%sop,0._8,in%wdir,.TRUE.)
CALL reporttime(0,0._8,in%reporttimefilename)
@@ -754,8 +765,8 @@ PROGRAM relax
in%x0,in%y0,in%lon0,in%lat0,in%zone,in%umult,in%wdir,oi)
END IF
#endif
+ WRITE (digit4,'(I4.4)') oi
#ifdef VTK
- WRITE (digit4,'(I4.4)') oi
IF (in%isoutputvtk .AND. in%isoutputstress) THEN
filename=trim(in%wdir)//"/sigma-"//digit4//".vtk"//char(0)
title="stress tensor field"//char(0)
@@ -770,6 +781,15 @@ PROGRAM relax
CALL exportvtk_rfaults_stress(in%sx1,in%sx2,in%sx3,in%dx1,in%dx2,in%dx3, &
in%nsop,in%sop,filename,convention=1,sig=sig)
#endif
+ ! total stress on predefined planes for gmt
+ filename=trim(in%wdir)//"/rfaults-sigma-"//digit4//".xy"
+ CALL exportgmt_rfaults_stress(in%sx1,in%sx2,in%sx3,in%dx1,in%dx2,in%dx3, &
+ in%nsop,in%sop,filename,sig=sig)
+ ! postseismic stress change on predefined planes for gm
+ filename=trim(in%wdir)//"/rfaults-dsigma-"//digit4//".xy"
+ CALL exportgmt_rfaults_stress(in%sx1,in%sx2,in%sx3,in%dx1,in%dx2,in%dx3, &
+ in%nsop,in%sop,filename,convention=1,sig=sig)
+ ! time series of stress in ASCII format
CALL exportcoulombstress(sig,in%sx1,in%sx2,in%sx3,in%dx1,in%dx2,in%dx3, &
in%nsop,in%sop,t,in%wdir,.FALSE.)
More information about the CIG-COMMITS
mailing list