[cig-commits] r22985 - seismo/3D/SPECFEM3D_GLOBE/trunk/src/auxiliaries
lefebvre at geodynamics.org
lefebvre at geodynamics.org
Thu Jan 2 08:00:41 PST 2014
Author: lefebvre
Date: 2014-01-02 08:00:41 -0800 (Thu, 02 Jan 2014)
New Revision: 22985
Modified:
seismo/3D/SPECFEM3D_GLOBE/trunk/src/auxiliaries/combine_vol_data.f90
seismo/3D/SPECFEM3D_GLOBE/trunk/src/auxiliaries/combine_vol_data_vtk.f90
Log:
combine_vol_data.f90 merged in combine_vol_data_vtk.f90
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/src/auxiliaries/combine_vol_data.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/src/auxiliaries/combine_vol_data.f90 2013-12-30 14:44:47 UTC (rev 22984)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/src/auxiliaries/combine_vol_data.f90 2014-01-02 16:00:41 UTC (rev 22985)
@@ -37,19 +37,28 @@
include "OUTPUT_FILES/values_from_mesher.h"
integer,parameter :: MAX_NUM_NODES = 2000
- integer iregion, ir, irs, ire, ires, pfd, efd
+ integer :: iregion, ir, irs, ire, ires
+!!! .mesh specific !!!!!!!!!!!
+ integer :: pfd, efd
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
character(len=256) :: sline, arg(7), filename, in_topo_dir, in_file_dir, outdir
character(len=256) :: prname_topo, prname_file, dimension_file
+!!! .mesh specific !!!!!!!!!!!
character(len=1038) :: command_name
- character(len=256) :: pt_mesh_file1, pt_mesh_file2, mesh_file, em_mesh_file, data_file, topo_file
+ character(len=256) :: pt_mesh_file1, pt_mesh_file2, mesh_file, em_mesh_file
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ character(len=256) :: data_file, topo_file
integer, dimension(MAX_NUM_NODES) :: node_list, nspec, nglob, npoint, nelement
integer iproc, num_node, i,j,k,ispec, ios, it, di, dj, dk
integer np, ne, njunk
+
real(kind=CUSTOM_REAL),dimension(NGLLX,NGLLY,NGLLZ,NSPEC_CRUST_MANTLE) :: data
real(kind=CUSTOM_REAL),dimension(NGLOB_CRUST_MANTLE) :: xstore, ystore, zstore
integer ibool(NGLLX,NGLLY,NGLLZ,NSPEC_CRUST_MANTLE)
+
integer num_ibool(NGLOB_CRUST_MANTLE)
logical mask_ibool(NGLOB_CRUST_MANTLE), HIGH_RESOLUTION_MESH
+
real x, y, z, dat
integer numpoin, iglob, n1, n2, n3, n4, n5, n6, n7, n8
integer iglob1, iglob2, iglob3, iglob4, iglob5, iglob6, iglob7, iglob8
@@ -175,6 +184,7 @@
do ir = irs, ire
print *, '----------- Region ', ir, '----------------'
+!!! .mesh specific !!!!!!!!!!!
! open paraview output mesh file
write(pt_mesh_file1,'(a,i1,a)') trim(outdir)//'/' // 'reg_',ir,'_'//trim(filename)//'_point1.mesh'
write(pt_mesh_file2,'(a,i1,a)') trim(outdir)//'/' // 'reg_',ir,'_'//trim(filename)//'_point2.mesh'
@@ -183,6 +193,7 @@
call open_file_fd(trim(pt_mesh_file1)//char(0),pfd)
call open_file_fd(trim(em_mesh_file)//char(0),efd)
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! figure out total number of points and elements for high-res mesh
@@ -206,6 +217,10 @@
read(27) nglob(it)
close(27)
+ ! check
+ if( nspec(it) > NSPEC_CRUST_MANTLE ) stop 'error file nspec too big, please check compilation'
+ if( nglob(it) > NGLOB_CRUST_MANTLE ) stop 'error file nglob too big, please check compilation'
+
if (HIGH_RESOLUTION_MESH) then
npoint(it) = nglob(it)
nelement(it) = nspec(it) * (NGLLX-1) * (NGLLY-1) * (NGLLZ-1)
@@ -231,6 +246,7 @@
iproc = node_list(it)
+ data(:,:,:,:) = -1.e10
print *, ' '
print *, 'Reading slice ', iproc
@@ -241,13 +257,16 @@
data_file = trim(prname_file) // trim(filename) // '.bin'
open(unit = 27,file = trim(data_file),status='old',action='read', iostat = ios,form ='unformatted')
if (ios /= 0) then
- print*,'error ',ios
- print*,'file:',trim(data_file)
- stop 'Error opening file'
+ print*,'error ',ios
+ print*,'file:',trim(data_file)
+ stop 'Error opening file'
endif
-
- data(:,:,:,:) = -1.e10
- read(27) data(:,:,:,1:nspec(it))
+ read(27,iostat=ios) data(:,:,:,1:nspec(it))
+ if( ios /= 0 ) then
+ print*,'read error ',ios
+ print*,'file:',trim(data_file)
+ stop 'error reading data'
+ endif
close(27)
print *,trim(data_file)
@@ -342,12 +361,12 @@
!dat = data(i,j,k,ispec)
dat = ibool_dat(iglob)
-
+!!! .mesh specific !!!!!!!!!!!
call write_real_fd(pfd,x)
call write_real_fd(pfd,y)
call write_real_fd(pfd,z)
call write_real_fd(pfd,dat)
-
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
mask_ibool(iglob) = .true.
num_ibool(iglob) = numpoin
endif
@@ -362,10 +381,12 @@
if( CORRECT_ELLIPTICITY ) call reverse_ellipticity(x,y,z,nspl,rspl,espl,espl2)
dat = data(i,j,k,ispec)
+!!! .mesh specific !!!!!!!!!!!
call write_real_fd(pfd,x)
call write_real_fd(pfd,y)
call write_real_fd(pfd,z)
call write_real_fd(pfd,dat)
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
mask_ibool(iglob) = .true.
num_ibool(iglob) = numpoin
endif
@@ -388,6 +409,7 @@
npoint(it) = numpoin
endif
+
! write elements file
numpoin = 0
do ispec = 1, nspec(it)
@@ -398,6 +420,7 @@
if( idoubling_inner_core(ispec) == IFLAG_IN_FICTITIOUS_CUBE) then
! connectivity must be given, otherwise element count would be wrong
! maps "fictitious" connectivity, element is all with iglob = 1
+!!! .mesh specific !!!!!!!!!!!
do k = 1, NGLLZ-1, dk
do j = 1, NGLLY-1, dj
do i = 1, NGLLX-1, di
@@ -412,16 +435,19 @@
enddo ! i
enddo ! j
enddo ! k
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! takes next element
cycle
endif
endif
! writes out element connectivity
- numpoin = numpoin + 1 ! counts elements
do k = 1, NGLLZ-1, dk
do j = 1, NGLLY-1, dj
do i = 1, NGLLX-1, di
+
+ numpoin = numpoin + 1 ! counts elements
+
iglob1 = ibool(i,j,k,ispec)
iglob2 = ibool(i+di,j,k,ispec)
iglob3 = ibool(i+di,j+dj,k,ispec)
@@ -430,6 +456,7 @@
iglob6 = ibool(i+di,j,k+dk,ispec)
iglob7 = ibool(i+di,j+dj,k+dk,ispec)
iglob8 = ibool(i,j+dj,k+dk,ispec)
+
n1 = num_ibool(iglob1)+np-1
n2 = num_ibool(iglob2)+np-1
n3 = num_ibool(iglob3)+np-1
@@ -438,6 +465,7 @@
n6 = num_ibool(iglob6)+np-1
n7 = num_ibool(iglob7)+np-1
n8 = num_ibool(iglob8)+np-1
+!!! .mesh specific !!!!!!!!!!!
call write_integer_fd(efd,n1)
call write_integer_fd(efd,n2)
call write_integer_fd(efd,n3)
@@ -446,6 +474,7 @@
call write_integer_fd(efd,n6)
call write_integer_fd(efd,n7)
call write_integer_fd(efd,n8)
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
enddo ! i
enddo ! j
enddo ! k
@@ -459,9 +488,11 @@
if (np /= sum(npoint(1:num_node))) stop 'Error: Number of total points are not consistent'
if (ne /= sum(nelement(1:num_node))) stop 'Error: Number of total elements are not consistent'
+ print *
print *, 'Total number of points: ', np
print *, 'Total number of elements: ', ne
-
+ print *
+!!! .mesh specific !!!!!!!!!!!
call close_file_fd(pfd)
call close_file_fd(efd)
@@ -475,7 +506,7 @@
print *, 'cat mesh files: '
print *, trim(command_name)
call system(trim(command_name))
-
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
enddo
print *, 'Done writing mesh files'
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/src/auxiliaries/combine_vol_data_vtk.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/src/auxiliaries/combine_vol_data_vtk.f90 2013-12-30 14:44:47 UTC (rev 22984)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/src/auxiliaries/combine_vol_data_vtk.f90 2014-01-02 16:00:41 UTC (rev 22985)
@@ -39,10 +39,9 @@
include "OUTPUT_FILES/values_from_mesher.h"
integer,parameter :: MAX_NUM_NODES = 2000
- integer iregion, ir, irs, ire, ires
+ integer :: iregion, ir, irs, ire, ires
character(len=256) :: sline, arg(7), filename, in_topo_dir, in_file_dir, outdir
character(len=256) :: prname_topo, prname_file, dimension_file
- character(len=256) :: mesh_file
character(len=256) :: data_file, topo_file
integer, dimension(MAX_NUM_NODES) :: node_list, nspec, nglob, npoint, nelement
integer iproc, num_node, i,j,k,ispec, ios, it, di, dj, dk
@@ -58,7 +57,6 @@
real x, y, z, dat
integer numpoin, iglob, n1, n2, n3, n4, n5, n6, n7, n8
integer iglob1, iglob2, iglob3, iglob4, iglob5, iglob6, iglob7, iglob8
- integer ier
! instead of taking the first value which appears for a global point, average the values
! if there are more than one gll points for a global point (points on element corners, edges, faces)
logical,parameter:: AVERAGE_GLOBALPOINTS = .false.
@@ -75,21 +73,33 @@
! puts point locations back into a perfectly spherical shape by removing the ellipticity factor;
! useful for plotting spherical cuts at certain depths
logical,parameter:: CORRECT_ELLIPTICITY = .false.
+
integer :: nspl
double precision :: rspl(NR),espl(NR),espl2(NR)
logical,parameter :: ONE_CRUST = .false. ! if you want to correct a model with one layer only in PREM crust
integer, dimension(NSPEC_INNER_CORE) :: idoubling_inner_core ! to get rid of fictitious elements in central cube
+!!! .mesh specific !!!!!!!!!!!
+#ifndef USE_VTK_INSTEAD_OF_MESH
+ integer :: pfd, efd
+ character(len=1038) :: command_name
+ character(len=256) :: pt_mesh_file1, pt_mesh_file2, mesh_file, em_mesh_file
+!!! .vtk specific !!!!!!!!!!!
+#else
+ character(len=256) :: mesh_file
+ integer ier
! global point data
real,dimension(:),allocatable :: total_dat
real,dimension(:,:),allocatable :: total_dat_xyz
integer,dimension(:,:),allocatable :: total_dat_con
+#endif
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! starts here--------------------------------------------------------------------------------------------------
do i = 1, 7
call get_command_argument(i,arg(i))
- if (i < 7 .and. trim(arg(i)) == '') then
+ if (i < 7 .and. len_trim(arg(i)) == 0) then
print *, ' '
print *, ' Usage: xcombine_vol_data slice_list filename input_topo_dir input_file_dir '
print *, ' output_dir high/low-resolution [region]'
@@ -107,7 +117,7 @@
stop 'This program needs that NSPEC_CRUST_MANTLE > NSPEC_OUTER_CORE and NSPEC_INNER_CORE'
! get region id
- if (trim(arg(7)) == '') then
+ if (len_trim(arg(7)) == 0) then
iregion = 0
else
read(arg(7),*) iregion
@@ -184,6 +194,18 @@
do ir = irs, ire
print *, '----------- Region ', ir, '----------------'
+!!! .mesh specific !!!!!!!!!!!
+#ifndef USE_VTK_INSTEAD_OF_MESH
+ ! open paraview output mesh file
+ write(pt_mesh_file1,'(a,i1,a)') trim(outdir)//'/' // 'reg_',ir,'_'//trim(filename)//'_point1.mesh'
+ write(pt_mesh_file2,'(a,i1,a)') trim(outdir)//'/' // 'reg_',ir,'_'//trim(filename)//'_point2.mesh'
+ write(mesh_file,'(a,i1,a)') trim(outdir)//'/' // 'reg_',ir,'_'//trim(filename)//'.mesh'
+ write(em_mesh_file,'(a,i1,a)') trim(outdir)//'/' // 'reg_',ir,'_'//trim(filename)//'_element.mesh'
+
+ call open_file_fd(trim(pt_mesh_file1)//char(0),pfd)
+ call open_file_fd(trim(em_mesh_file)//char(0),efd)
+#endif
+
! figure out total number of points and elements for high-res mesh
do it = 1, num_node
@@ -214,10 +236,16 @@
npoint(it) = nglob(it)
nelement(it) = nspec(it) * (NGLLX-1) * (NGLLY-1) * (NGLLZ-1)
else if( ires == 0 ) then
+!!! .vtk specific !!!!!!!!!!!
+#ifdef USE_VTK_INSTEAD_OF_MESH
npoint(it) = nglob(it)
+#endif
nelement(it) = nspec(it)
else if (ires == 2 ) then
+!!! .vtk specific !!!!!!!!!!!
+#ifdef USE_VTK_INSTEAD_OF_MESH
npoint(it) = nglob(it)
+#endif
nelement(it) = nspec(it) * (NGLLX-1) * (NGLLY-1) * (NGLLZ-1) / 8
endif
@@ -225,9 +253,11 @@
print *, 'nspec(it) = ', nspec(1:num_node)
print *, 'nglob(it) = ', nglob(1:num_node)
+ print *, 'nelement(it) = ', nelement(1:num_node)
- !call write_integer_fd(efd,sum(nelement(1:num_node)))
-
+#ifndef USE_VTK_INSTEAD_OF_MESH
+ call write_integer_fd(efd,sum(nelement(1:num_node)))
+#else
! VTK
print *
print *,'vtk inital total points: ',sum(npoint(1:num_node))
@@ -244,6 +274,7 @@
allocate(total_dat_con(8,sum(nelement(1:num_node))),stat=ier)
if( ier /= 0 ) stop 'error allocating total_dat_con array'
total_dat_con(:,:) = 0
+#endif
np = 0
ne = 0
@@ -340,7 +371,6 @@
num_ibool(:) = 0
numpoin = 0
-
! write point file
do ispec=1,nspec(it)
! checks if element counts
@@ -370,18 +400,18 @@
!dat = data(i,j,k,ispec)
dat = ibool_dat(iglob)
-
- !call write_real_fd(pfd,x)
- !call write_real_fd(pfd,y)
- !call write_real_fd(pfd,z)
- !call write_real_fd(pfd,dat)
-
+#ifndef USE_VTK_INSTEAD_OF_MESH
+ call write_real_fd(pfd,x)
+ call write_real_fd(pfd,y)
+ call write_real_fd(pfd,z)
+ call write_real_fd(pfd,dat)
+#else
! VTK
total_dat(np+numpoin) = dat
total_dat_xyz(1,np+numpoin) = x
total_dat_xyz(2,np+numpoin) = y
total_dat_xyz(3,np+numpoin) = z
-
+#endif
mask_ibool(iglob) = .true.
num_ibool(iglob) = numpoin
endif
@@ -397,17 +427,18 @@
dat = data(i,j,k,ispec)
- !call write_real_fd(pfd,x)
- !call write_real_fd(pfd,y)
- !call write_real_fd(pfd,z)
- !call write_real_fd(pfd,dat)
-
+#ifndef USE_VTK_INSTEAD_OF_MESH
+ call write_real_fd(pfd,x)
+ call write_real_fd(pfd,y)
+ call write_real_fd(pfd,z)
+ call write_real_fd(pfd,dat)
+#else
! VTK
total_dat(np+numpoin) = dat
total_dat_xyz(1,np+numpoin) = x
total_dat_xyz(2,np+numpoin) = y
total_dat_xyz(3,np+numpoin) = z
-
+#endif
mask_ibool(iglob) = .true.
num_ibool(iglob) = numpoin
endif
@@ -417,7 +448,6 @@
enddo ! k
enddo !ispec
-
! no way to check the number of points for low-res
if (HIGH_RESOLUTION_MESH ) then
if( ir==3 ) then
@@ -441,20 +471,22 @@
if( idoubling_inner_core(ispec) == IFLAG_IN_FICTITIOUS_CUBE) then
! connectivity must be given, otherwise element count would be wrong
! maps "fictitious" connectivity, element is all with iglob = 1
- !do k = 1, NGLLZ-1, dk
- ! do j = 1, NGLLY-1, dj
- ! do i = 1, NGLLX-1, di
- !call write_integer_fd(efd,1)
- !call write_integer_fd(efd,1)
- !call write_integer_fd(efd,1)
- !call write_integer_fd(efd,1)
- !call write_integer_fd(efd,1)
- !call write_integer_fd(efd,1)
- !call write_integer_fd(efd,1)
- !call write_integer_fd(efd,1)
- ! enddo ! i
- ! enddo ! j
- !enddo ! k
+#ifndef USE_VTK_INSTEAD_OF_MESH
+ do k = 1, NGLLZ-1, dk
+ do j = 1, NGLLY-1, dj
+ do i = 1, NGLLX-1, di
+ call write_integer_fd(efd,1)
+ call write_integer_fd(efd,1)
+ call write_integer_fd(efd,1)
+ call write_integer_fd(efd,1)
+ call write_integer_fd(efd,1)
+ call write_integer_fd(efd,1)
+ call write_integer_fd(efd,1)
+ call write_integer_fd(efd,1)
+ enddo ! i
+ enddo ! j
+ enddo ! k
+#endif
! takes next element
cycle
endif
@@ -485,15 +517,16 @@
n7 = num_ibool(iglob7)+np-1
n8 = num_ibool(iglob8)+np-1
- !call write_integer_fd(efd,n1)
- !call write_integer_fd(efd,n2)
- !call write_integer_fd(efd,n3)
- !call write_integer_fd(efd,n4)
- !call write_integer_fd(efd,n5)
- !call write_integer_fd(efd,n6)
- !call write_integer_fd(efd,n7)
- !call write_integer_fd(efd,n8)
-
+#ifndef USE_VTK_INSTEAD_OF_MESH
+ call write_integer_fd(efd,n1)
+ call write_integer_fd(efd,n2)
+ call write_integer_fd(efd,n3)
+ call write_integer_fd(efd,n4)
+ call write_integer_fd(efd,n5)
+ call write_integer_fd(efd,n6)
+ call write_integer_fd(efd,n7)
+ call write_integer_fd(efd,n8)
+#else
! VTK
! note: indices for vtk start at 0
total_dat_con(1,numpoin + ne) = n1
@@ -504,7 +537,7 @@
total_dat_con(6,numpoin + ne) = n6
total_dat_con(7,numpoin + ne) = n7
total_dat_con(8,numpoin + ne) = n8
-
+#endif
enddo ! i
enddo ! j
enddo ! k
@@ -523,6 +556,7 @@
print *, 'Total number of elements: ', ne
print *
+#ifdef USE_VTK_INSTEAD_OF_MESH
! VTK
! opens unstructured grid file
write(mesh_file,'(a,i1,a)') trim(outdir)//'/' // 'reg_',ir,'_'//trim(filename)//'.vtk'
@@ -548,21 +582,6 @@
total_dat_con(5,i),total_dat_con(6,i),total_dat_con(7,i),total_dat_con(8,i)
enddo
write(IOVTK,*) ""
-
- !call close_file_fd(pfd)
- !call close_file_fd(efd)
-
- ! add the critical piece: total number of points
- !call open_file_fd(trim(pt_mesh_file2)//char(0),pfd)
- !call write_integer_fd(pfd,np)
- !call close_file_fd(pfd)
-
- !command_name='cat '//trim(pt_mesh_file2)//' '//trim(pt_mesh_file1)//' '//trim(em_mesh_file)//' > '//trim(mesh_file)
- !print *, ' '
- !print *, 'cat mesh files: '
- !print *, trim(command_name)
- !call system(trim(command_name))
-
! VTK
! type: hexahedrons
write(IOVTK,'(a,i12)') "CELL_TYPES ",ne
@@ -584,6 +603,21 @@
print *,'written: ',trim(mesh_file)
print *
+#else
+ call close_file_fd(pfd)
+ call close_file_fd(efd)
+
+ ! add the critical piece: total number of points
+ call open_file_fd(trim(pt_mesh_file2)//char(0),pfd)
+ call write_integer_fd(pfd,np)
+ call close_file_fd(pfd)
+
+ command_name='cat '//trim(pt_mesh_file2)//' '//trim(pt_mesh_file1)//' '//trim(em_mesh_file)//' > '//trim(mesh_file)
+ print *, ' '
+ print *, 'cat mesh files: '
+ print *, trim(command_name)
+ call system(trim(command_name))
+#endif
enddo
print *, 'Done writing mesh files'
More information about the CIG-COMMITS
mailing list