adds vtk write routine

Author: danielpeter

src/shared/write_VTK_file.f90

@@ -1474,10 +1474,108 @@ end subroutine write_VTK_movie_data_elemental
 !  close(IOUT_VTK)
 !
 !  end subroutine write_VTK_movie_data_binary
+
 !
+!-------------------------------------------------------------------------------------------------
 !
+
+  subroutine write_VTK_wavefield(nspec,nglob,xstore,ystore,zstore,ibool, &
+                                 field,prname_file)
+
+  use constants, only: CUSTOM_REAL,MAX_STRING_LEN,IOUT_VTK,NDIM,NGLLX,NGLLY,NGLLZ
+
+  implicit none
+
+  integer,intent(in) :: nspec,nglob
+
+  ! global coordinates
+  integer, dimension(NGLLX,NGLLY,NGLLZ,nspec),intent(in) :: ibool
+  real(kind=CUSTOM_REAL), dimension(nglob),intent(in) :: xstore,ystore,zstore
+
+  ! GLL data values array
+  real(kind=CUSTOM_REAL), dimension(NDIM,nglob),intent(in) :: field
+
+  ! file name
+  character(len=MAX_STRING_LEN),intent(in) :: prname_file
+
+  ! local parameters
+  integer :: ispec,i,ier
+
+  open(IOUT_VTK,file=trim(prname_file)//'.vtk',status='unknown',iostat=ier)
+  if (ier /= 0 ) then
+    print *, 'Error opening VTK output file: ',trim(prname_file)
+    stop 'Error opening VTK output file'
+  endif
+  write(IOUT_VTK,'(a)') '# vtk DataFile Version 3.1'
+  write(IOUT_VTK,'(a)') 'material model VTK file'
+  write(IOUT_VTK,'(a)') 'ASCII'
+  write(IOUT_VTK,'(a)') 'DATASET UNSTRUCTURED_GRID'
+  write(IOUT_VTK, '(a,i12,a)') 'POINTS ', nglob, ' float'
+  do i = 1,nglob
+    write(IOUT_VTK,'(3e18.6)') real(xstore(i),kind=4),real(ystore(i),kind=4),real(zstore(i),kind=4)
+  enddo
+  write(IOUT_VTK,*) ""
+
+  ! note: indices for vtk start at 0
+  write(IOUT_VTK,'(a,i12,i12)') "CELLS ",nspec,nspec*9
+  do ispec = 1,nspec
+    write(IOUT_VTK,'(9i12)') 8, &
+                             ibool(1,1,1,ispec)-1,ibool(NGLLX,1,1,ispec)-1, &
+                             ibool(NGLLX,NGLLY,1,ispec)-1,ibool(1,NGLLY,1,ispec)-1, &
+                             ibool(1,1,NGLLZ,ispec)-1,ibool(NGLLX,1,NGLLZ,ispec)-1, &
+                             ibool(NGLLX,NGLLY,NGLLZ,ispec)-1,ibool(1,NGLLY,NGLLZ,ispec)-1
+  enddo
+  write(IOUT_VTK,*) ""
+
+  ! type: hexahedrons
+  write(IOUT_VTK,'(a,i12)') "CELL_TYPES ",nspec
+  write(IOUT_VTK,'(6i12)') (12,ispec=1,nspec)
+  write(IOUT_VTK,*) ""
+
+  write(IOUT_VTK,'(a,i12)') "POINT_DATA ",nglob
+
+  ! single wavefield components
+  if (.true.) then
+    write(IOUT_VTK,'(a)') "SCALARS field_x float"
+    write(IOUT_VTK,'(a)') "LOOKUP_TABLE default"
+    do i = 1,nglob
+      write(IOUT_VTK,*) real(field(1,i),kind=4)
+    enddo
+    write(IOUT_VTK,*) ""
+
+    write(IOUT_VTK,'(a)') "SCALARS field_y float"
+    write(IOUT_VTK,'(a)') "LOOKUP_TABLE default"
+    do i = 1,nglob
+      write(IOUT_VTK,*) real(field(2,i),kind=4)
+    enddo
+    write(IOUT_VTK,*) ""
+
+    write(IOUT_VTK,'(a)') "SCALARS field_z float"
+    write(IOUT_VTK,'(a)') "LOOKUP_TABLE default"
+    do i = 1,nglob
+      write(IOUT_VTK,*) real(field(3,i),kind=4)
+    enddo
+    write(IOUT_VTK,*) ""
+  endif
+
+  ! vector wavefield
+  if (.true.) then
+    write(IOUT_VTK,'(a)') "VECTORS field_vector float"
+    do i = 1,nglob
+      write(IOUT_VTK,*) real(field(1,i),kind=4),real(field(2,i),kind=4),real(field(3,i),kind=4)
+    enddo
+    write(IOUT_VTK,*) ""
+  endif
+
+  close(IOUT_VTK)
+
+  end subroutine write_VTK_wavefield
+
 !-------------------------------------------------------------------------------------------------
 !
+! VTU binary formats
+!
+!-------------------------------------------------------------------------------------------------
 
   subroutine write_VTU_movie_data(ne,np,total_dat_xyz,total_dat_con,total_dat,mesh_file,var_name)