[cig-commits] commit: added Coulomb stress calculation, added comments in example runs, and export of coulomb stress in vtk
Mercurial
hg at geodynamics.org
Tue Sep 20 12:13:19 PDT 2011
changeset: 20:4b0d45863b7f
user: Sylvain Barbot <sylbar.vainbot at gmail.com>
date: Tue Aug 09 19:19:51 2011 -0700
files: examples/run1.sh examples/run2.sh examples/run3.sh export.f90 input.f90 makefile_imkl relax.f90 types.f90
description:
added Coulomb stress calculation, added comments in example runs, and export of coulomb stress in vtk
diff -r 4a5123723424 -r 4b0d45863b7f examples/run1.sh
--- a/examples/run1.sh Tue May 31 10:18:15 2011 -0700
+++ b/examples/run1.sh Tue Aug 09 19:19:51 2011 -0700
@@ -66,7 +66,7 @@ if [ ! -e $WDIR ]; then
mkdir $WDIR
fi
-time ../relax $* <<EOF | tee output1/in.param
+time ../relax $* <<EOF | tee $WDIR/in.param
# use '#' character to include comments in your input file
# grid size (sx1,sx2,sx3)
256 256 256
@@ -91,6 +91,8 @@ 20 -1 1
# number of observation planes
0
# number of observation points
+0
+# number of Coulomb patches
0
# number of prestress interfaces
0
diff -r 4a5123723424 -r 4b0d45863b7f examples/run2.sh
--- a/examples/run2.sh Tue May 31 10:18:15 2011 -0700
+++ b/examples/run2.sh Tue Aug 09 19:19:51 2011 -0700
@@ -23,27 +23,45 @@ if [ ! -e $WDIR ]; then
mkdir $WDIR
fi
-time ../relax --no-proj-output $* <<EOF | tee output2/in.param
+time ../relax --no-proj-output $* <<EOF | tee $WDIR/in.param
# grid size (sx1,sx2,sx3)
256 256 256
+# sampling size defines the grid spacing in units of distance
+# smoothing is a roll-off parameter 0 <= beta <= 0.5 that tapers
+# the slip distribution on faults beta=0.2 and beta=0.3 are good values
+# nyquist is the minimum size of fault patches, relative to the sampling size
# sampling size, smoothing & nyquist (dx1,dx2,dx3,beta,nq)
0.05 0.05 0.05 0.2 2
# origin position & rotation
0 0 0
# observation depth (displacements and stress) (stress is only exported in GRD)
0 0.5
-# output directory
+# output directory (all output files are stored in the following directory)
$WDIR
+# the elastic parameters (Lambda and mu) are the Lame parameter and the shear modulus
+# gamma=(1-nu)*rho*g/mu is a wavelength relevant to buoyancy.
# elastic parameters and gamma = (1-nu) rho g / mu = 8.33e-7 /m = 8.33e-4 /km
1 1 8.33e-4
+# integration time refers to the length in units of time of the simulation.
# integration time (t1)
-100 -2 1
+100 -1 1
+# observation planes are arbitrary planes where the inelastic strain rate
+# is sampled for output.
# number of observation planes
0
+# observation points are GPS points, for example, where time series of
+# displacement are output.
# number of observation points
0
+# stress observation segments are patches where Coulomb and other
+# stress components are evaluated.
+# number of stress observation segments
+0
+# prestress interfaces corresponds to depths where the prestress
+# tensor changes (default is no prestress).
# number of prestress interfaces
0
+# viscous interfaces are depths where the viscous properties change.
# number of linear viscous interfaces
0
# number of powerlaw viscous interfaces
@@ -51,23 +69,27 @@ 2
# no depth gammadot0 power cohesion
1 3.0 5e3 3.0 0.0
2 8.0 5e3 3.0 0.0
+# ductile zones corresponds to volumes where viscous properties change.
# number of power-law ductile zones
0
-# number of friction faults
+# depths where friction properties change.
+# number of friction interfaces
0
# number of interseismic loading strike-slip and opening
0
0
# number of coseismic events
1
-# number of shear dislocations
+# number of shear dislocations (strike-slip and dip-slip faults)
1
# index slip x1 x2 x3 length width strike dip rake
1 1 -1 0 0 2 1 0 90 0
# number of tensile cracks
0
+# volumetric deformation (Mogi source)
# number of dilatation sources
0
+# surface loads can be used to model loading of lakes or dams
# number of surface loads
0
EOF
diff -r 4a5123723424 -r 4b0d45863b7f examples/run3.sh
--- a/examples/run3.sh Tue May 31 10:18:15 2011 -0700
+++ b/examples/run3.sh Tue Aug 09 19:19:51 2011 -0700
@@ -34,7 +34,7 @@ if [ ! -e $WDIR ]; then
mkdir $WDIR
fi
-time ../relax --no-proj-output --no-stress-output $* <<EOF | tee output3/in.param
+time ../relax --no-proj-output --no-stress-output $* <<EOF | tee $WDIR/in.param
# grid size (sx1,sx2,sx3)
256 256 256
# sampling size, smoothing & nyquist (dx1,dx2,dx3,beta,nq)
@@ -59,13 +59,15 @@ 0
0
# number of observation points
0
+# number of Coulomb patches
+0
# number of prestress interfaces
0
# number of linear viscous interfaces
0
# number of powerlaw viscous interfaces
0
-# number of friction faults
+# number of friction interfaces
1
# no depth gamma0 (a-b)sig friction cohesion
1 1 1e3 1e3 0.6 0
diff -r 4a5123723424 -r 4b0d45863b7f export.f90
--- a/export.f90 Tue May 31 10:18:15 2011 -0700
+++ b/export.f90 Tue Aug 09 19:19:51 2011 -0700
@@ -962,9 +962,9 @@ END SUBROUTINE exportcreep
!!
!! \author sylvain barbot 06/24/09 - original form
!------------------------------------------------------------------
- SUBROUTINE exportvtk_grid(sx1,sx2,sx3,dx1,dx2,dx3,origx,origy,cgfilename)
+ SUBROUTINE exportvtk_grid(sx1,sx2,sx3,dx1,dx2,dx3,cgfilename)
INTEGER, INTENT(IN) :: sx1,sx2,sx3
- REAL*8, INTENT(IN) :: dx1,dx2,dx3,origx,origy
+ REAL*8, INTENT(IN) :: dx1,dx2,dx3
CHARACTER(80), INTENT(IN) :: cgfilename
INTEGER :: iostatus
@@ -1220,6 +1220,241 @@ END SUBROUTINE exportcreep
CLOSE(15)
END SUBROUTINE exportvtk_rfaults
+
+ !------------------------------------------------------------------
+ !> subroutine ExportVTK_RFaults_Stress_Init
+ !! 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 exportvtk_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
+ REAL*8, INTENT(IN) :: dx1,dx2,dx3
+ TYPE(SEGMENT_STRUCT), INTENT(INOUT), DIMENSION(nsop) :: sop
+
+ INTEGER :: k,i1,i2,i3
+ REAL*8 :: x1,x2,x3
+ ! local value of stress
+ TYPE(TENSOR) :: lsig
+
+ DO k=1,nsop
+ ! fault center position
+ x1=sop(k)%x
+ x2=sop(k)%y
+ x3=sop(k)%z
+
+ CALL shiftedindex(x1,x2,x3,sx1,sx2,sx3,dx1,dx2,dx3,i1,i2,i3)
+ lsig=sig(i1,i2,i3)
+
+ sop(k)%sig0%s11=lsig%s11
+ sop(k)%sig0%s12=lsig%s12
+ sop(k)%sig0%s13=lsig%s13
+ sop(k)%sig0%s22=lsig%s22
+ sop(k)%sig0%s23=lsig%s23
+ sop(k)%sig0%s33=lsig%s33
+
+ END DO
+
+ END SUBROUTINE exportvtk_rfaults_stress_init
+
+ !------------------------------------------------------------------
+ !> subroutine ExportVTK_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 exportvtk_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
+
+ ! 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
+
+ WRITE (15,'("<?xml version=",a,"1.0",a,"?>")') q,q
+ WRITE (15,'("<VTKFile type=",a,"PolyData",a," version=",a,"0.1",a,">")') q,q,q,q
+ WRITE (15,'(" <PolyData>")')
+
+ 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,'(" <Piece NumberOfPoints=",a,"4",a," NumberOfPolys=",a,"1",a,">")'),q,q,q,q
+ WRITE (15,'(" <Points>")')
+ WRITE (15,'(" <DataArray type=",a,"Float32",a, &
+ " Name=",a,"Fault Patch",a, &
+ " NumberOfComponents=",a,"3",a, &
+ " 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
+ WRITE (15,'(" </DataArray>")')
+
+ WRITE (15,'(" </Points>")')
+ WRITE (15,'(" <Polys>")')
+ WRITE (15,'(" <DataArray type=",a,"Int32",a, &
+ " Name=",a,"connectivity",a, &
+ " format=",a,"ascii",a, &
+ " RangeMin=",a,"0",a, &
+ " RangeMax=",a,"3",a,">")'), q,q,q,q,q,q,q,q,q,q
+ WRITE (15,'("0 1 2 3")')
+ WRITE (15,'(" </DataArray>")')
+ WRITE (15,'(" <DataArray type=",a,"Int32",a, &
+ " Name=",a,"offsets",a, &
+ " format=",a,"ascii",a, &
+ " RangeMin=",a,"4",a, &
+ " RangeMax=",a,"4",a,">")'), q,q,q,q,q,q,q,q,q,q
+ WRITE (15,'(" 4")')
+ WRITE (15,'(" </DataArray>")')
+ WRITE (15,'(" </Polys>")')
+
+ WRITE (15,'(" <CellData Normals=",a,"stress",a,">")'), q,q
+
+ WRITE (15,'(" <DataArray type=",a,"Float32",a, &
+ " Name=",a,"stress tensor",a, &
+ " NumberOfComponents=",a,"6",a, &
+ " format=",a,"ascii",a,">")'), q,q,q,q,q,q,q,q
+ WRITE (15,'(6ES11.2)'), lsig%s11,lsig%s12,lsig%s13,lsig%s22,lsig%s23,lsig%s33
+ WRITE (15,'(" </DataArray>")')
+
+ WRITE (15,'(" <DataArray type=",a,"Float32",a, &
+ " Name=",a,"shear stress",a, &
+ " NumberOfComponents=",a,"1",a, &
+ " format=",a,"ascii",a,">")'), q,q,q,q,q,q,q,q
+ WRITE (15,'(ES11.2)'), taus
+ WRITE (15,'(" </DataArray>")')
+
+ WRITE (15,'(" <DataArray type=",a,"Float32",a, &
+ " Name=",a,"normal stress",a, &
+ " NumberOfComponents=",a,"1",a, &
+ " format=",a,"ascii",a,">")'), q,q,q,q,q,q,q,q
+ WRITE (15,'(ES11.2)'), taun
+ WRITE (15,'(" </DataArray>")')
+
+ WRITE (15,'(" </CellData>")')
+
+ WRITE (15,'(" </Piece>")')
+
+ END DO
+
+ WRITE (15,'(" </PolyData>")')
+ WRITE (15,'("</VTKFile>")')
+
+ CLOSE(15)
+
+ END SUBROUTINE exportvtk_rfaults_stress
!------------------------------------------------------------------
!> subroutine ExportVTK_Rectangle
@@ -1588,7 +1823,7 @@ END SUBROUTINE exportcreep
REAL*8, INTENT(IN) :: dx1,dx2,dx3
CHARACTER(80), INTENT(IN) :: vcfilename,title,name
- INTEGER :: iostatus,idum,i1,i2,i3
+ INTEGER :: iostatus,idum,i1,i2
CHARACTER :: q,end_rec
CHARACTER(LEN=180) :: s_buffer
REAL*8 :: x1,x2,x3
@@ -1657,7 +1892,7 @@ END SUBROUTINE exportcreep
REAL*8, INTENT(IN) :: dx1,dx2,dx3
CHARACTER(80), INTENT(IN) :: vcfilename,title,name
- INTEGER :: iostatus,idum,i1,i2,i3
+ INTEGER :: iostatus,idum,i1,i2
CHARACTER :: q,end_rec
CHARACTER(LEN=180) :: s_buffer
REAL*8 :: x1,x2,x3
diff -r 4a5123723424 -r 4b0d45863b7f input.f90
--- a/input.f90 Tue May 31 10:18:15 2011 -0700
+++ b/input.f90 Tue Aug 09 19:19:51 2011 -0700
@@ -48,10 +48,10 @@ CONTAINS
CHARACTER :: ch
CHARACTER(180) :: dataline
- CHARACTER(80) :: inputfilename
- CHARACTER(80) :: rffilename,cgfilename
+ CHARACTER(80) :: rffilename,filename
#ifdef VTK
CHARACTER(3) :: digit
+ CHARACTER(4) :: digit4
#endif
INTEGER :: iunit
!$ INTEGER :: omp_get_num_procs,omp_get_max_threads
@@ -237,7 +237,6 @@ CONTAINS
in%reporttimefilename=trim(in%wdir)//"/time.txt"
in%reportfilename=trim(in%wdir)//"/report.txt"
- inputfilename=trim(in%wdir)//"/relax.inp"
#ifdef TXT
PRINT '(" ",a," (report: ",a,")")', trim(in%wdir),trim(in%reportfilename)
#else
@@ -256,9 +255,8 @@ CONTAINS
! end test
#ifdef VTK
- cgfilename=trim(in%wdir)//"/cgrid.vtp"
- CALL exportvtk_grid(in%sx1,in%sx2,in%sx3,in%dx1,in%dx2,in%dx3, &
- in%x0,in%y0,cgfilename)
+ filename=trim(in%wdir)//"/cgrid.vtp"
+ CALL exportvtk_grid(in%sx1,in%sx2,in%sx3,in%dx1,in%dx2,in%dx3,filename)
#endif
PRINT '(a)', "lambda, mu, gamma (gamma = (1 - nu) rho g / mu)"
@@ -352,6 +350,61 @@ CONTAINS
STOP 1
END IF
END DO
+
+ END IF
+
+ ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
+ ! C O U L O M B O B S E R V A T I O N S E G M E N T S
+ ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
+ PRINT '(a)', "number of stress observation segments"
+ CALL getdata(iunit,dataline)
+ READ (dataline,*) in%nsop
+ PRINT '(I5)', in%nsop
+ IF (in%nsop .gt. 0) THEN
+ ALLOCATE(in%sop(in%nsop),STAT=iostatus)
+ IF (iostatus>0) STOP "could not allocate the segment list"
+ PRINT 2000
+ PRINT '(a)',"no. xs ys zs length width strike dip friction"
+ PRINT 2000
+ DO k=1,in%nsop
+ CALL getdata(iunit,dataline)
+ READ (dataline,*) i, &
+ in%sop(k)%x,in%sop(k)%y,in%sop(k)%z, &
+ in%sop(k)%length,in%sop(k)%width, &
+ in%sop(k)%strike,in%sop(k)%dip,in%sop(k)%friction
+ in%sop(k)%sig0=TENSOR(0.d0,0.d0,0.d0,0.d0,0.d0,0.d0)
+
+ PRINT '(I4.4,3ES9.2E1,2ES8.2E1,f7.1,f6.1,f7.1)',i, &
+ in%sop(k)%x,in%sop(k)%y,in%sop(k)%z, &
+ in%sop(k)%length,in%sop(k)%width, &
+ in%sop(k)%strike,in%sop(k)%dip, &
+ in%sop(k)%friction
+
+ IF (i .ne. k) THEN
+ WRITE_DEBUG_INFO
+ WRITE (0,'("invalid segment definition ")')
+ WRITE (0,'("error in input file: source index misfit")')
+ STOP 1
+ END IF
+ IF (MAX(in%sop(k)%length,in%sop(k)%width) .LE. 0._8) THEN
+ WRITE_DEBUG_INFO
+ WRITE (0,'("error in input file: length and width must be positive.")')
+ STOP 1
+ END IF
+
+ ! comply to Wang's convention
+ CALL wangconvention(dummy, &
+ in%sop(k)%x,in%sop(k)%y,in%sop(k)%z, &
+ in%sop(k)%length,in%sop(k)%width, &
+ in%sop(k)%strike,in%sop(k)%dip, &
+ dummy, &
+ in%x0,in%y0,in%rot)
+ END DO
+
+ ! export patches to vtk/vtp
+ filename=trim(in%wdir)//"/rfaults-dsigma-0000.vtp"
+ CALL exportvtk_rfaults_stress(in%sx1,in%sx2,in%sx3,in%dx1,in%dx2,in%dx3, &
+ in%nsop,in%sop,filename,convention=1)
END IF
@@ -682,9 +735,9 @@ CONTAINS
#ifdef VTK
! export the afterslip segment in VTK format
- WRITE (digit,'(I3.3)') k
+ WRITE (digit4,'(I4.4)') k
- rffilename=trim(in%wdir)//"/aplane-"//digit//".vtp"
+ rffilename=trim(in%wdir)//"/aplane-"//digit4//".vtp"
CALL exportvtk_rectangle(in%n(k)%x,in%n(k)%y,in%n(k)%z, &
in%n(k)%length,in%n(k)%width, &
in%n(k)%strike,in%n(k)%dip,rffilename)
diff -r 4a5123723424 -r 4b0d45863b7f makefile_imkl
--- a/makefile_imkl Tue May 31 10:18:15 2011 -0700
+++ b/makefile_imkl Tue Aug 09 19:19:51 2011 -0700
@@ -14,7 +14,7 @@ INCPATH=-I/sw/include -I/opt/intel/Compi
LIBS=-lproj -lgmt -lnetcdf -lmkl_intel_lp64 -lmkl_intel_thread -lmkl_core -openmp -lpthread
-CC=cc
+CC=cc -em64
F77=ifort
FC=ifort -openmp
diff -r 4a5123723424 -r 4b0d45863b7f relax.f90
--- a/relax.f90 Tue May 31 10:18:15 2011 -0700
+++ b/relax.f90 Tue Aug 09 19:19:51 2011 -0700
@@ -210,8 +210,9 @@ PROGRAM relax
REAL*8 :: maxwell(3)
TYPE(SIMULATION_STRUC) :: in
#ifdef VTK
- CHARACTER(80) :: vcfilename,title,name
+ CHARACTER(80) :: filename,title,name
CHARACTER(3) :: digit
+ CHARACTER(4) :: digit4
#endif
REAL*8 :: t,Dt,tm
@@ -321,7 +322,8 @@ PROGRAM relax
! test the presence of dislocations for coseismic calculation
IF ((in%events(e)%nt .NE. 0) .OR. &
(in%events(e)%ns .NE. 0) .OR. &
- (in%events(e)%nm .NE. 0)) THEN
+ (in%events(e)%nm .NE. 0) .OR. &
+ (in%events(e)%nl .NE. 0)) THEN
! apply the 3d elastic transfer function
CALL greenfunctioncowling(v1,v2,v3,t1,t2,t3, &
@@ -359,14 +361,14 @@ PROGRAM relax
#endif
#ifdef VTK
IF (in%isoutputvtk) THEN
- !vcfilename=trim(in%wdir)//"/disp-000.vtr"
- !CALL exportvtk_vectors(u1,u2,u3,in%sx1,in%sx2,in%sx3/4,in%dx1,in%dx2,in%dx3,8,8,8,vcfilename)
- vcfilename=trim(in%wdir)//"/disp-000.vtk"
+ !filename=trim(in%wdir)//"/disp-000.vtr"
+ !CALL exportvtk_vectors(u1,u2,u3,in%sx1,in%sx2,in%sx3/4,in%dx1,in%dx2,in%dx3,8,8,8,filename)
+ filename=trim(in%wdir)//"/disp-000.vtk"
title="coseismic displacement vector field"
name="displacement"
CALL exportvtk_vectors_legacy(u1,u2,u3,in%sx1,in%sx2,in%sx3/4,in%dx1,in%dx2,in%dx3, &
- 4,4,8,vcfilename,title,name)
- !CALL exportvtk_vectors_slice(u1,u2,u3,in%sx1,in%sx2,in%sx3/2,in%dx1,in%dx2,in%dx3,in%oz,8,8,vcfilename)
+ 4,4,8,filename,title,name)
+ !CALL exportvtk_vectors_slice(u1,u2,u3,in%sx1,in%sx2,in%sx3/2,in%dx1,in%dx2,in%dx3,in%oz,8,8,filename)
END IF
#endif
IF (ALLOCATED(in%ptsname)) THEN
@@ -413,12 +415,24 @@ PROGRAM relax
#endif
#ifdef VTK
IF (in%isoutputvtk) THEN
- WRITE (digit,'(I3.3)') oi-1
- vcfilename=trim(in%wdir)//"/sigma-"//digit//".vtk"
+ WRITE (digit4,'(I4.4)') oi-1
+ filename=trim(in%wdir)//"/sigma-"//digit4//".vtk"
title="stress tensor field"
name="stress"
CALL exportvtk_tensors_legacy(sig,in%sx1,in%sx2,in%sx3/2,in%dx1,in%dx2,in%dx3, &
- 4,4,8,vcfilename,title,name)
+ 4,4,8,filename,title,name)
+ END IF
+ 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)
+ IF (0.EQ.(oi-1)) THEN
+ ! 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)
+ ELSE
+ 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,convention=1,sig=sig)
END IF
#endif
END IF
@@ -467,11 +481,11 @@ PROGRAM relax
#ifdef VTK
IF (in%isoutputvtk) THEN
WRITE (digit,'(I3.3)') oi-1
- vcfilename=trim(in%wdir)//"/power-"//digit//".vtk"
+ filename=trim(in%wdir)//"/power-"//digit//".vtk"
title="stress rate tensor field"
name="power"
CALL exportvtk_tensors_legacy(moment,in%sx1,in%sx2,in%sx3/2,in%dx1,in%dx2,in%dx3, &
- 4,4,8,vcfilename,title,name)
+ 4,4,8,filename,title,name)
END IF
#endif
@@ -575,13 +589,13 @@ PROGRAM relax
#ifdef VTK_EQBF
IF (in%isoutputvtk) THEN
WRITE (digit,'(I3.3)') oi
- !vcfilename=trim(in%wdir)//"/eqbf-"//digit//".vtr"
- !CALL exportvtk_vectors(v1,v2,v3,in%sx1,in%sx2,in%sx3/4,in%dx1,in%dx2,in%dx3,8,8,8,vcfilename)
- vcfilename=trim(in%wdir)//"/eqbf-"//digit//".vtk"
+ !filename=trim(in%wdir)//"/eqbf-"//digit//".vtr"
+ !CALL exportvtk_vectors(v1,v2,v3,in%sx1,in%sx2,in%sx3/4,in%dx1,in%dx2,in%dx3,8,8,8,filename)
+ filename=trim(in%wdir)//"/eqbf-"//digit//".vtk"
title="instantaneous equivalent body-force rate vector field"
name="body-force-rate"
CALL exportvtk_vectors_legacy(v1,v2,v3,in%sx1,in%sx2,in%sx3/4,in%dx1,in%dx2,in%dx3, &
- 4,4,8,vcfilename,title,name)
+ 4,4,8,filename,title,name)
END IF
#endif
#ifdef GRD_EQBF
@@ -680,24 +694,24 @@ PROGRAM relax
IF (in%isoutputvtk) THEN
WRITE (digit,'(I3.3)') oi
! export total displacement in VTK XML format
- !vcfilename=trim(in%wdir)//"/disp-"//digit//".vtr"
- !CALL exportvtk_vectors(u1,u2,u3,in%sx1,in%sx2,in%sx3/4,in%dx1,in%dx2,in%dx3,8,8,8,vcfilename)
- vcfilename=trim(in%wdir)//"/disp-"//digit//".vtk"
+ !filename=trim(in%wdir)//"/disp-"//digit//".vtr"
+ !CALL exportvtk_vectors(u1,u2,u3,in%sx1,in%sx2,in%sx3/4,in%dx1,in%dx2,in%dx3,8,8,8,filename)
+ filename=trim(in%wdir)//"/disp-"//digit//".vtk"
title="cumulative displacement vector field"
name="displacement"
CALL exportvtk_vectors_legacy(u1,u2,u3,in%sx1,in%sx2,in%sx3/4,in%dx1,in%dx2,in%dx3, &
- 4,4,8,vcfilename,title,name)
- !CALL exportvtk_vectors_slice(u1,u2,u3,in%sx1,in%sx2,in%sx3/2,in%dx1,in%dx2,in%dx3,in%oz,8,8,vcfilename)
+ 4,4,8,filename,title,name)
+ !CALL exportvtk_vectors_slice(u1,u2,u3,in%sx1,in%sx2,in%sx3/2,in%dx1,in%dx2,in%dx3,in%oz,8,8,filename)
! export instantaneous velocity in VTK XML format
- !vcfilename=trim(in%wdir)//"/vel-"//digit//".vtr"
- !CALL exportvtk_vectors(v1,v2,v3,in%sx1,in%sx2,in%sx3/4,in%dx1,in%dx2,in%dx3,8,8,8,vcfilename)
- vcfilename=trim(in%wdir)//"/vel-"//digit//".vtk"
+ !filename=trim(in%wdir)//"/vel-"//digit//".vtr"
+ !CALL exportvtk_vectors(v1,v2,v3,in%sx1,in%sx2,in%sx3/4,in%dx1,in%dx2,in%dx3,8,8,8,filename)
+ filename=trim(in%wdir)//"/vel-"//digit//".vtk"
title="instantaneous velocity vector field"
name="velocity"
CALL exportvtk_vectors_legacy(v1,v2,v3,in%sx1,in%sx2,in%sx3/2,in%dx1,in%dx2,in%dx3, &
- 8,8,16,vcfilename,title,name)
- !CALL exportvtk_vectors_slice(v1,v2,v3,in%sx1,in%sx2,in%sx3/2,in%dx1,in%dx2,in%dx3,in%oz,8,8,vcfilename)
+ 8,8,16,filename,title,name)
+ !CALL exportvtk_vectors_slice(v1,v2,v3,in%sx1,in%sx2,in%sx3/2,in%dx1,in%dx2,in%dx3,in%oz,8,8,filename)
END IF
#endif
@@ -728,6 +742,7 @@ 100 CONTINUE
IF (ALLOCATED(in%opts)) DEALLOCATE(in%opts)
IF (ALLOCATED(in%ptsname)) DEALLOCATE(in%ptsname)
IF (ALLOCATED(in%op)) DEALLOCATE(in%op)
+ IF (ALLOCATED(in%sop)) DEALLOCATE(in%sop)
IF (ALLOCATED(in%n)) DEALLOCATE(in%n)
IF (ALLOCATED(in%stressstruc)) DEALLOCATE(in%stressstruc)
IF (ALLOCATED(in%stresslayer)) DEALLOCATE(in%stresslayer)
diff -r 4a5123723424 -r 4b0d45863b7f types.f90
--- a/types.f90 Tue May 31 10:18:15 2011 -0700
+++ b/types.f90 Tue Aug 09 19:19:51 2011 -0700
@@ -57,6 +57,12 @@ MODULE types
TYPE(TENSOR) :: t
END TYPE TENSOR_LAYER_STRUCT
+ TYPE SEGMENT_STRUCT
+ SEQUENCE
+ REAL*8 :: x,y,z,width,length,strike,dip,friction
+ TYPE(TENSOR) :: sig0
+ END TYPE SEGMENT_STRUCT
+
TYPE SLIPPATCH_STRUCT
SEQUENCE
REAL*8 :: x1,x2,x3,lx,lz,slip,ss,ds
@@ -112,6 +118,12 @@ MODULE types
! observation planes
TYPE(PLANE_STRUCT), DIMENSION(:), ALLOCATABLE :: op
+
+ ! number of stress observation planes
+ INTEGER :: nsop
+
+ ! stress observation planes
+ TYPE(SEGMENT_STRUCT), DIMENSION(:), ALLOCATABLE :: sop
! number of observation points
INTEGER :: npts
More information about the CIG-COMMITS
mailing list