[cig-commits] [commit] master: add traction component output for GMT. (99c95fb)
cig_noreply at geodynamics.org
cig_noreply at geodynamics.org
Wed Apr 23 21:26:44 PDT 2014
Repository : ssh://geoshell/relax
On branch : master
Link : https://github.com/geodynamics/relax/compare/24d6b5c21fe72c085a3e8dfb56f1a445d0f85490...99c95fb11bdd3448888f349ba2a0c5dd712b1e71
>---------------------------------------------------------------
commit 99c95fb11bdd3448888f349ba2a0c5dd712b1e71
Author: Sylvain Barbot <sbarbot at ntu.edu.sg>
Date: Thu Apr 24 12:26:20 2014 +0800
add traction component output for GMT.
>---------------------------------------------------------------
99c95fb11bdd3448888f349ba2a0c5dd712b1e71
src/export.f90 | 110 +++++++++++++++++++++++++++++++++++++++++++++++++++++++++
src/relax.f90 | 3 ++
2 files changed, 113 insertions(+)
diff --git a/src/export.f90 b/src/export.f90
index 42ee987..d43d6a8 100644
--- a/src/export.f90
+++ b/src/export.f90
@@ -1692,6 +1692,116 @@ END SUBROUTINE exportcreep_vtk
END SUBROUTINE export_rfaults_stress_init
!------------------------------------------------------------------
+ !> subroutine ExportGMT_RFaults_Traction
+ !! creates a .dat file containing fault center coordinates and
+ !! horizontal traction components.
+ !!
+ !! \author sylvain barbot 04/23/14 - original form
+ !------------------------------------------------------------------
+ SUBROUTINE exportgmt_rfaults_traction(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(256), 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,taun
+ ! segment normal vector, strike direction, dip direction
+ REAL*8, DIMENSION(3) :: n,s,d
+ ! local value of stress
+ TYPE(TENSOR) :: lsig
+ ! 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
+
+ ! 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
+
+ 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
+
+ WRITE (15,'(5ES11.2)') x2,x1,ts(2),ts(1)
+
+ END DO
+
+ CLOSE(15)
+
+ END SUBROUTINE exportgmt_rfaults_traction
+
+ !------------------------------------------------------------------
!> 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
diff --git a/src/relax.f90 b/src/relax.f90
index 1edd3d5..227f74c 100644
--- a/src/relax.f90
+++ b/src/relax.f90
@@ -441,6 +441,9 @@ PROGRAM relax
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)
+ filename=trim(in%wdir)//"/rfaults-traction-"//digit4//".dat"
+ CALL exportgmt_rfaults_traction(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, &
More information about the CIG-COMMITS
mailing list