[cig-commits] commit: fixes output index of stress field, export stress tensor field in VTK (binary).
Mercurial
hg at geodynamics.org
Tue Sep 20 12:13:18 PDT 2011
changeset: 19:4a5123723424
user: Sylvain Barbot <sylbar.vainbot at gmail.com>
date: Tue May 31 10:18:15 2011 -0700
files: .doxygen export.f90 relax.f90
description:
fixes output index of stress field, export stress tensor field in VTK (binary).
diff -r 88603e5baf07 -r 4a5123723424 .doxygen
--- a/.doxygen Mon May 23 21:48:51 2011 -0700
+++ b/.doxygen Tue May 31 10:18:15 2011 -0700
@@ -867,7 +867,7 @@ ENUM_VALUES_PER_LINE = 4
# releases of Doxygen, the values YES and NO are equivalent to FRAME and NONE
# respectively.
-GENERATE_TREEVIEW = NONE
+GENERATE_TREEVIEW = YES
# If the treeview is enabled (see GENERATE_TREEVIEW) then this tag can be
# used to set the initial width (in pixels) of the frame in which the tree
diff -r 88603e5baf07 -r 4a5123723424 export.f90
--- a/export.f90 Mon May 23 21:48:51 2011 -0700
+++ b/export.f90 Tue May 31 10:18:15 2011 -0700
@@ -1644,6 +1644,77 @@ END SUBROUTINE exportcreep
END SUBROUTINE exportvtk_vectors_legacy
!------------------------------------------------------------------
+ !> subroutine ExportVTK_tensors_Legacy
+ !! creates a .vtk file in the VTK Legacy binary format with
+ !! structured points containing a vector field.
+ !!
+ !! \author sylvain barbot 05/23/11 - original form
+ !------------------------------------------------------------------
+ SUBROUTINE exportvtk_tensors_legacy(sig,sx1,sx2,sx3,dx1,dx2,dx3, &
+ j1,j2,j3,vcfilename,title,name)
+ INTEGER, INTENT(IN) :: sx1,sx2,sx3,j1,j2,j3
+ TYPE(TENSOR), INTENT(IN), DIMENSION(sx1,sx2,sx3) :: sig
+ REAL*8, INTENT(IN) :: dx1,dx2,dx3
+ CHARACTER(80), INTENT(IN) :: vcfilename,title,name
+
+ INTEGER :: iostatus,idum,i1,i2,i3
+ CHARACTER :: q,end_rec
+ CHARACTER(LEN=180) :: s_buffer
+ REAL*8 :: x1,x2,x3
+
+ ! double-quote character
+ q=char(34)
+ ! end-character for binary-record finalize
+ end_rec=char(10)
+
+ ! open file for mixed binary/ascii writing in VTK legacy format
+ OPEN(UNIT=15,FILE=vcfilename,form='UNFORMATTED',ACCESS='SEQUENTIAL', &
+ ACTION='WRITE',CONVERT='BIG_ENDIAN',RECORDTYPE='STREAM', &
+ BUFFERED='YES',IOSTAT=iostatus)
+ IF (iostatus>0) THEN
+ WRITE_DEBUG_INFO
+ PRINT '(a)', vcfilename
+ STOP "could not open file in mixed binary/ascii for export in legacy format"
+ END IF
+
+ ! writing header of file (INIT)
+ WRITE(UNIT=15,IOSTAT=iostatus) '# vtk DataFile Version 3.0'//end_rec
+ WRITE(UNIT=15,IOSTAT=iostatus) TRIM(title)//end_rec
+ WRITE(UNIT=15,IOSTAT=iostatus) 'BINARY'//end_rec
+ WRITE(UNIT=15,IOSTAT=iostatus) 'DATASET '//'STRUCTURED_POINTS'//end_rec
+
+ ! structured points grid (GEO)
+ WRITE(s_buffer,FMT='(A,3I12)',IOSTAT=iostatus) 'DIMENSIONS ',sx1/j1,sx2/j2,sx3/j3
+ WRITE(UNIT=15,IOSTAT=iostatus) TRIM(s_buffer)//end_rec
+ WRITE(s_buffer,FMT='(A,3E14.6E2)',IOSTAT=iostatus) 'ORIGIN ',-dx1*(sx1/2),-dx2*(sx2/2),0.d0
+ WRITE(UNIT=15,IOSTAT=iostatus) TRIM(s_buffer)//end_rec
+ WRITE(s_buffer,FMT='(A,3E14.6E2)',IOSTAT=iostatus) 'SPACING ',dx1*j1,dx2*j2,dx3*j3
+ WRITE(UNIT=15,IOSTAT=iostatus) TRIM(s_buffer)//end_rec
+
+ ! data header for this grid (DAT)
+ WRITE(s_buffer,FMT='(A,I12)',IOSTAT=iostatus) 'POINT_DATA ', (sx1/j1)*(sx2/j2)*(sx3/j3)
+ WRITE(UNIT=15,IOSTAT=iostatus) TRIM(s_buffer)//end_rec
+
+ ! data array (VAR)
+ WRITE(UNIT=15,IOSTAT=iostatus) 'TENSORS '//TRIM(name)//' float'//end_rec
+ DO x3=0,sx3-1,j3
+ DO x2=-sx2/2,sx2/2-1,j2
+ DO x1=-sx1/2,sx1/2-1,j1
+ CALL shiftedindex(x1,x2,1._8,sx1,sx2,sx3,1._8,1._8,1._8,i1,i2,idum)
+ WRITE (UNIT=15,IOSTAT=iostatus) sig(i1,i2,x3+1)%s11,sig(i1,i2,x3+1)%s12,sig(i1,i2,x3+1)%s13, &
+ sig(i1,i2,x3+1)%s12,sig(i1,i2,x3+1)%s22,sig(i1,i2,x3+1)%s23, &
+ sig(i1,i2,x3+1)%s13,sig(i1,i2,x3+1)%s23,sig(i1,i2,x3+1)%s33
+ END DO
+ END DO
+ END DO
+ WRITE(UNIT=15,IOSTAT=iostatus) end_rec
+
+ ! close binary file (END)
+ CLOSE(15)
+
+ END SUBROUTINE exportvtk_tensors_legacy
+
+ !------------------------------------------------------------------
!> subroutine ExportVTK_Vectors_Slice
!! creates a .vtr file (in the VTK Rectilinear XML format)
!! containing a vector field.
diff -r 88603e5baf07 -r 4a5123723424 relax.f90
--- a/relax.f90 Mon May 23 21:48:51 2011 -0700
+++ b/relax.f90 Tue May 31 10:18:15 2011 -0700
@@ -158,22 +158,22 @@
!!
!! MODIFICATIONS:
!! \author Sylvain Barbot
- !! (07-06-07) - original form
+ !! (07-06-07) - original form <br>
!! (08-28-08) - FFTW/SGI_FFT support, FIR derivatives,
!! Runge-Kutta integration, tensile cracks,
- !! GMT output, comments in input file
+ !! GMT output, comments in input file <br>
!! (10-24-08) - interseismic loading, postseismic signal
- !! output in separate files
- !! (12-08-09) - slip distribution smoothing
+ !! output in separate files <br>
+ !! (12-08-09) - slip distribution smoothing <br>
!! (05-05-10) - lateral variations in viscous properties
- !! Intel MKL implementation of the FFT
+ !! Intel MKL implementation of the FFT <br>
!! (06-04-10) - output in geographic coordinates
- !! and output components of stress tensor
+ !! and output components of stress tensor <br>
!! (07-19-10) - includes surface tractions initial condition
- !! output geometry in VTK format for Paraview
+ !! output geometry in VTK format for Paraview <br>
!! (02-28-11) - add constraints on the broad direction of
!! afterslip, export faults to GMT xy format
- !! and allow scaling of computed time steps.
+ !! and allow scaling of computed time steps. <br>
!! (04-26-11) - include command-line arguments
!!
!! \todo
@@ -182,6 +182,8 @@
!! - evaluate Green function, stress and body forces in GPU
!! - create a ./configure ./install framework
!! - input a list of planes to compute Coulomb stress
+ !! - export a .xy list of observation points with names
+ !! - fix output index
!------------------------------------------------------------------------
PROGRAM relax
@@ -400,13 +402,23 @@ PROGRAM relax
#ifdef GRD
IF (in%isoutputgrd .AND. in%isoutputstress) THEN
CALL exportstressgrd(sig,in%sx1,in%sx2,in%sx3/2,in%dx1,in%dx2,in%dx3, &
- in%ozs,in%x0,in%y0,in%wdir,i-1)
+ in%ozs,in%x0,in%y0,in%wdir,oi-1)
END IF
#endif
#ifdef PROJ
IF (in%isoutputproj .AND. in%isoutputstress) THEN
CALL exportstressproj(sig,in%sx1,in%sx2,in%sx3/2,in%dx1,in%dx2,in%dx3,in%ozs, &
- in%x0,in%y0,in%lon0,in%lat0,in%zone,in%umult,in%wdir,i-1)
+ in%x0,in%y0,in%lon0,in%lat0,in%zone,in%umult,in%wdir,oi-1)
+ END IF
+#endif
+#ifdef VTK
+ IF (in%isoutputvtk) THEN
+ WRITE (digit,'(I3.3)') oi-1
+ vcfilename=trim(in%wdir)//"/sigma-"//digit//".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)
END IF
#endif
END IF
@@ -451,6 +463,17 @@ PROGRAM relax
END DO
mech(3)=1
END IF
+
+#ifdef VTK
+ IF (in%isoutputvtk) THEN
+ WRITE (digit,'(I3.3)') oi-1
+ vcfilename=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)
+ END IF
+#endif
! identify the required time step
tm=1._8/(REAL(mech(1))/maxwell(1)+ &
More information about the CIG-COMMITS
mailing list