[cig-commits] r20599 - in seismo/3D/SPECFEM3D/trunk/src: generate_databases meshfem3D shared specfem3D
dkomati1 at geodynamics.org
dkomati1 at geodynamics.org
Mon Aug 20 05:57:46 PDT 2012
Author: dkomati1
Date: 2012-08-20 05:57:46 -0700 (Mon, 20 Aug 2012)
New Revision: 20599
Modified:
seismo/3D/SPECFEM3D/trunk/src/generate_databases/create_regions_mesh.f90
seismo/3D/SPECFEM3D/trunk/src/meshfem3D/check_mesh_quality.f90
seismo/3D/SPECFEM3D/trunk/src/meshfem3D/create_regions_mesh.f90
seismo/3D/SPECFEM3D/trunk/src/meshfem3D/define_subregions.f90
seismo/3D/SPECFEM3D/trunk/src/meshfem3D/define_subregions_heuristic.f90
seismo/3D/SPECFEM3D/trunk/src/meshfem3D/meshfem3D.f90
seismo/3D/SPECFEM3D/trunk/src/meshfem3D/store_boundaries.f90
seismo/3D/SPECFEM3D/trunk/src/meshfem3D/store_coords.f90
seismo/3D/SPECFEM3D/trunk/src/shared/check_mesh_resolution.f90
seismo/3D/SPECFEM3D/trunk/src/shared/constants.h.in
seismo/3D/SPECFEM3D/trunk/src/shared/create_movie_shakemap_AVS_DX_GMT.f90
seismo/3D/SPECFEM3D/trunk/src/specfem3D/setup_movie_meshes.f90
seismo/3D/SPECFEM3D/trunk/src/specfem3D/specfem3D_par.f90
seismo/3D/SPECFEM3D/trunk/src/specfem3D/write_movie_output.f90
Log:
started to add support for HEX27 elements; not finished yet
Modified: seismo/3D/SPECFEM3D/trunk/src/generate_databases/create_regions_mesh.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/generate_databases/create_regions_mesh.f90 2012-08-20 12:15:45 UTC (rev 20598)
+++ seismo/3D/SPECFEM3D/trunk/src/generate_databases/create_regions_mesh.f90 2012-08-20 12:57:46 UTC (rev 20599)
@@ -638,7 +638,6 @@
zstore(:,:,:,:) = 0.d0
do ispec = 1, nspec
- !call get_xyzelm(xelm, yelm, zelm, ispec, elmnts_ext_mesh, nodes_coords_ext_mesh, nspec, nnodes_ext_mesh)
do ia = 1,NGNOD
xelm(ia) = nodes_coords_ext_mesh(1,elmnts_ext_mesh(ia,ispec))
yelm(ia) = nodes_coords_ext_mesh(2,elmnts_ext_mesh(ia,ispec))
Modified: seismo/3D/SPECFEM3D/trunk/src/meshfem3D/check_mesh_quality.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/meshfem3D/check_mesh_quality.f90 2012-08-20 12:15:45 UTC (rev 20598)
+++ seismo/3D/SPECFEM3D/trunk/src/meshfem3D/check_mesh_quality.f90 2012-08-20 12:57:46 UTC (rev 20599)
@@ -24,7 +24,7 @@
!
!=====================================================================
-! read a 2D or 3D CUBIT mesh file and display statistics about mesh quality;
+! read a 3D CUBIT mesh file and display statistics about mesh quality;
! and create an OpenDX file showing a given range of elements or a single element
@@ -50,7 +50,7 @@
integer :: NSPEC
double precision, dimension(NGLOB) :: x,y,z
- integer, dimension(NGNOD,NSPEC) :: ibool
+ integer, dimension(NGNOD_EIGHT_CORNERS,NSPEC) :: ibool
logical :: CREATE_VTK_FILES
character(len=256) prname
@@ -231,14 +231,9 @@
write(IMAIN,*) '*** computed using VP_MAX = ',VP_MAX
write(IMAIN,*) '***'
- ! max stability CFL value is different in 2D and in 3D
- if(NGNOD == 8) then
- max_CFL_stability_limit = 0.48d0
- else if(NGNOD == 4) then
- max_CFL_stability_limit = 0.68d0
- else
- stop 'NGNOD must be 4 or 8'
- endif
+ ! max stability CFL value
+!! DK DK we could probably increase this, now that the Stacey conditions have been fixed
+ max_CFL_stability_limit = 0.55d0 !! DK DK increased this 0.48d0
if(stability_max_MPI >= max_CFL_stability_limit) then
write(IMAIN,*) '*********************************************'
@@ -396,9 +391,9 @@
double precision, dimension(NGLOB) :: x,y,z
- integer, dimension(NGNOD,NSPEC) :: ibool
+ integer, dimension(NGNOD_EIGHT_CORNERS,NSPEC) :: ibool
- double precision, dimension(NGNOD) :: xelm,yelm,zelm
+ double precision, dimension(NGNOD_EIGHT_CORNERS) :: xelm,yelm,zelm
double precision vectorA_x,vectorA_y,vectorA_z
double precision vectorB_x,vectorB_y,vectorB_z
@@ -417,7 +412,7 @@
integer faces_topo(6,6)
! store the corners of this element for the skewness routine
- do i = 1,NGNOD
+ do i = 1,NGNOD_EIGHT_CORNERS
xelm(i) = x(ibool(i,ispec))
yelm(i) = y(ibool(i,ispec))
zelm(i) = z(ibool(i,ispec))
Modified: seismo/3D/SPECFEM3D/trunk/src/meshfem3D/create_regions_mesh.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/meshfem3D/create_regions_mesh.f90 2012-08-20 12:15:45 UTC (rev 20598)
+++ seismo/3D/SPECFEM3D/trunk/src/meshfem3D/create_regions_mesh.f90 2012-08-20 12:57:46 UTC (rev 20599)
@@ -101,13 +101,13 @@
integer ner_layer(nblayers)
! topology of the elements
- integer iaddx(NGNOD)
- integer iaddy(NGNOD)
- integer iaddz(NGNOD)
+ integer iaddx(NGNOD_EIGHT_CORNERS)
+ integer iaddy(NGNOD_EIGHT_CORNERS)
+ integer iaddz(NGNOD_EIGHT_CORNERS)
- double precision xelm(NGNOD)
- double precision yelm(NGNOD)
- double precision zelm(NGNOD)
+ double precision xelm(NGNOD_EIGHT_CORNERS)
+ double precision yelm(NGNOD_EIGHT_CORNERS)
+ double precision zelm(NGNOD_EIGHT_CORNERS)
! boundary locator
logical, dimension(:,:), allocatable :: iboun
@@ -266,8 +266,8 @@
! Regular subregion case
- ! loop over the NGNOD nodes
- do ia=1,NGNOD
+ ! loop over the NGNOD_EIGHT_CORNERS nodes
+ do ia=1,NGNOD_EIGHT_CORNERS
! define topological coordinates of this mesh point
ioffset_x = ix+iax*iaddx(ia)
ioffset_y = iy+iay*iaddy(ia)
Modified: seismo/3D/SPECFEM3D/trunk/src/meshfem3D/define_subregions.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/meshfem3D/define_subregions.f90 2012-08-20 12:15:45 UTC (rev 20598)
+++ seismo/3D/SPECFEM3D/trunk/src/meshfem3D/define_subregions.f90 2012-08-20 12:57:46 UTC (rev 20599)
@@ -42,9 +42,9 @@
integer num_material
! topology of the elements
- integer iaddx(NGNOD)
- integer iaddy(NGNOD)
- integer iaddz(NGNOD)
+ integer iaddx(NGNOD_EIGHT_CORNERS)
+ integer iaddy(NGNOD_EIGHT_CORNERS)
+ integer iaddz(NGNOD_EIGHT_CORNERS)
integer ner_layer(nblayers)
! definition of the different regions of the model in the mesh (nx,ny,nz)
@@ -110,9 +110,9 @@
integer ndoublings
! topology of the elements
- integer iaddx(NGNOD)
- integer iaddy(NGNOD)
- integer iaddz(NGNOD)
+ integer iaddx(NGNOD_EIGHT_CORNERS)
+ integer iaddy(NGNOD_EIGHT_CORNERS)
+ integer iaddz(NGNOD_EIGHT_CORNERS)
integer ner_layer(nblayers)
integer ner_doublings(2)
Modified: seismo/3D/SPECFEM3D/trunk/src/meshfem3D/define_subregions_heuristic.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/meshfem3D/define_subregions_heuristic.f90 2012-08-20 12:15:45 UTC (rev 20598)
+++ seismo/3D/SPECFEM3D/trunk/src/meshfem3D/define_subregions_heuristic.f90 2012-08-20 12:57:46 UTC (rev 20599)
@@ -45,9 +45,9 @@
integer NER_BOTTOM_MOHO,NER_MOHO_16,NER_16_BASEMENT,NER_BASEMENT_SEDIM
! topology of the elements
- integer iaddx(NGNOD)
- integer iaddy(NGNOD)
- integer iaddz(NGNOD)
+ integer iaddx(NGNOD_EIGHT_CORNERS)
+ integer iaddy(NGNOD_EIGHT_CORNERS)
+ integer iaddz(NGNOD_EIGHT_CORNERS)
! **************
Modified: seismo/3D/SPECFEM3D/trunk/src/meshfem3D/meshfem3D.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/meshfem3D/meshfem3D.f90 2012-08-20 12:15:45 UTC (rev 20598)
+++ seismo/3D/SPECFEM3D/trunk/src/meshfem3D/meshfem3D.f90 2012-08-20 12:57:46 UTC (rev 20599)
@@ -516,8 +516,9 @@
if(CUSTOM_REAL /= SIZE_REAL .and. CUSTOM_REAL /= SIZE_DOUBLE) call exit_MPI(myrank,'wrong size of CUSTOM_REAL for reals')
! checks number of nodes for hexahedra and quadrilaterals
- if(NGNOD /= 8) call exit_MPI(myrank,'number of control nodes must be 8')
- if(NGNOD2D /= 4) call exit_MPI(myrank,'elements with 8 points should have NGNOD2D = 4')
+ ! curvature (i.e., HEX27 elements) is not handled by our internal mesher, for that use Gmsh (CUBIT does not handle it either)
+ if(NGNOD /= NGNOD_EIGHT_CORNERS) call exit_MPI(myrank,'number of control nodes must be 8 in our internal mesher')
+ if(NGNOD2D /= NGNOD2D_FOUR_CORNERS) call exit_MPI(myrank,'elements should have NGNOD2D = 4 in our internal mesher')
! for the number of standard linear solids for attenuation
if(N_SLS /= 3) call exit_MPI(myrank,'number of SLS must be 3')
Modified: seismo/3D/SPECFEM3D/trunk/src/meshfem3D/store_boundaries.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/meshfem3D/store_boundaries.f90 2012-08-20 12:15:45 UTC (rev 20598)
+++ seismo/3D/SPECFEM3D/trunk/src/meshfem3D/store_boundaries.f90 2012-08-20 12:57:46 UTC (rev 20599)
@@ -50,9 +50,8 @@
! counters to keep track of number of elements on each of the boundaries
integer ispecb1,ispecb2,ispecb3,ispecb4,ispecb5,ispecb6
- ! check that the parameter file is correct
- if(NGNOD /= 8) call exit_MPI(myrank,'elements should have 8 control nodes')
- if(NGNOD2D /= 4) call exit_MPI(myrank,'surface elements should have 4 control nodes')
+ ! check that the parameter is correct
+ if(NGNOD2D_FOUR_CORNERS /= 4) call exit_MPI(myrank,'surface elements should have 4 control nodes in our internal mesher')
! initializes
ispecb1 = 0
@@ -171,10 +170,10 @@
real(kind=CUSTOM_REAL) normal_bottom(NDIM,NGLLX,NGLLY,NSPEC2D_BOTTOM)
real(kind=CUSTOM_REAL) normal_top(NDIM,NGLLX,NGLLY,NSPEC2D_TOP)
- double precision dershape2D_x(NDIM2D,NGNOD2D,NGLLY,NGLLZ)
- double precision dershape2D_y(NDIM2D,NGNOD2D,NGLLX,NGLLZ)
- double precision dershape2D_bottom(NDIM2D,NGNOD2D,NGLLX,NGLLY)
- double precision dershape2D_top(NDIM2D,NGNOD2D,NGLLX,NGLLY)
+ double precision dershape2D_x(NDIM2D,NGNOD2D_FOUR_CORNERS,NGLLY,NGLLZ)
+ double precision dershape2D_y(NDIM2D,NGNOD2D_FOUR_CORNERS,NGLLX,NGLLZ)
+ double precision dershape2D_bottom(NDIM2D,NGNOD2D_FOUR_CORNERS,NGLLX,NGLLY)
+ double precision dershape2D_top(NDIM2D,NGNOD2D_FOUR_CORNERS,NGLLX,NGLLY)
! global element numbering
integer ispec
@@ -182,11 +181,10 @@
! counters to keep track of number of elements on each of the boundaries
integer ispecb1,ispecb2,ispecb3,ispecb4,ispecb5,ispecb6
- double precision xelm(NGNOD2D),yelm(NGNOD2D),zelm(NGNOD2D)
+ double precision xelm(NGNOD2D_FOUR_CORNERS),yelm(NGNOD2D_FOUR_CORNERS),zelm(NGNOD2D_FOUR_CORNERS)
-! check that the parameter file is correct
- if(NGNOD /= 8) call exit_MPI(myrank,'elements should have 8 control nodes')
- if(NGNOD2D /= 4) call exit_MPI(myrank,'surface elements should have 4 control nodes')
+ ! check that the parameter is correct
+ if(NGNOD2D_FOUR_CORNERS /= 4) call exit_MPI(myrank,'surface elements should have 4 control nodes in our internal mesher')
ispecb1 = 0
ispecb2 = 0
@@ -380,8 +378,8 @@
integer ispecb,NGLLA,NGLLB,NSPEC2DMAX_AB,myrank
- double precision xelm(NGNOD2D),yelm(NGNOD2D),zelm(NGNOD2D)
- double precision dershape2D(NDIM2D,NGNOD2D,NGLLA,NGLLB)
+ double precision xelm(NGNOD2D_FOUR_CORNERS),yelm(NGNOD2D_FOUR_CORNERS),zelm(NGNOD2D_FOUR_CORNERS)
+ double precision dershape2D(NDIM2D,NGNOD2D_FOUR_CORNERS,NGLLA,NGLLB)
real(kind=CUSTOM_REAL) jacobian2D(NGLLA,NGLLB,NSPEC2DMAX_AB)
real(kind=CUSTOM_REAL) normal(3,NGLLA,NGLLB,NSPEC2DMAX_AB)
@@ -399,7 +397,7 @@
yeta=ZERO
zxi=ZERO
zeta=ZERO
- do ia=1,NGNOD2D
+ do ia=1,NGNOD2D_FOUR_CORNERS
xxi=xxi+dershape2D(1,ia,i,j)*xelm(ia)
xeta=xeta+dershape2D(2,ia,i,j)*xelm(ia)
yxi=yxi+dershape2D(1,ia,i,j)*yelm(ia)
Modified: seismo/3D/SPECFEM3D/trunk/src/meshfem3D/store_coords.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/meshfem3D/store_coords.f90 2012-08-20 12:15:45 UTC (rev 20598)
+++ seismo/3D/SPECFEM3D/trunk/src/meshfem3D/store_coords.f90 2012-08-20 12:57:46 UTC (rev 20599)
@@ -32,7 +32,7 @@
integer ispec,nspec
- double precision, dimension(NGNOD) :: xelm,yelm,zelm
+ double precision, dimension(NGNOD_EIGHT_CORNERS) :: xelm,yelm,zelm
double precision, dimension(NGLLX,NGLLY,NGLLZ,nspec) :: xstore,ystore,zstore
Modified: seismo/3D/SPECFEM3D/trunk/src/shared/check_mesh_resolution.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/shared/check_mesh_resolution.f90 2012-08-20 12:15:45 UTC (rev 20598)
+++ seismo/3D/SPECFEM3D/trunk/src/shared/check_mesh_resolution.f90 2012-08-20 12:57:46 UTC (rev 20599)
@@ -960,8 +960,8 @@
! corners indices of reference cube faces
! shapes of arrays below
- integer,dimension(2),parameter :: corner_shape = (/3,NGNOD/)
- integer,dimension(3,NGNOD),parameter :: corner_ijk = &
+ integer,dimension(2),parameter :: corner_shape = (/3,NGNOD_EIGHT_CORNERS/)
+ integer,dimension(3,NGNOD_EIGHT_CORNERS),parameter :: corner_ijk = &
reshape((/ 1,1,1, 1,NGLLY,1, 1,NGLLY,NGLLZ, 1,1,NGLLZ, &
NGLLX,1,1, NGLLX,NGLLY,1, NGLLX,NGLLY,NGLLZ, NGLLX,1,NGLLZ /),corner_shape)
@@ -970,7 +970,7 @@
elemsize_max = -HUGEVAL
! loops over corners
- do icorner=1,NGNOD
+ do icorner=1,NGNOD_EIGHT_CORNERS
i = corner_ijk(1,icorner)
j = corner_ijk(2,icorner)
k = corner_ijk(3,icorner)
@@ -981,7 +981,7 @@
z0 = zstore(iglob_a)
! loops over all other corners
- do jcorner = icorner+1,NGNOD
+ do jcorner = icorner+1,NGNOD_EIGHT_CORNERS
i = corner_ijk(1,jcorner)
j = corner_ijk(2,jcorner)
k = corner_ijk(3,jcorner)
Modified: seismo/3D/SPECFEM3D/trunk/src/shared/constants.h.in
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/shared/constants.h.in 2012-08-20 12:15:45 UTC (rev 20598)
+++ seismo/3D/SPECFEM3D/trunk/src/shared/constants.h.in 2012-08-20 12:57:46 UTC (rev 20599)
@@ -45,6 +45,14 @@
!----------- parameters that can be changed by the user -----------
+! number of nodes for 2D and 3D shape functions for hexahedra
+! we use either 8-node mesh elements (bricks) or 27-node elements.
+! If you use our internal mesher, the only option is 8-node bricks (27-node elements are not supported)
+! CUBIT does not support HEX27 elements either (it can generate them, but they are flat, i.e. identical to HEX8).
+! To generate HEX27 elements with curvature properly taken into account, you can use Gmsh http://geuz.org/gmsh/
+ integer, parameter :: NGNOD = 8, NGNOD2D = 4
+! integer, parameter :: NGNOD = 27, NGNOD2D = 9
+
! set to .false. if running on a Beowulf-type machine with local disks
! set to .true. if running on a shared-memory machine with common file system
! if running on a Beowulf, also modify name of nodes in filter_machine_file.f90
@@ -288,13 +296,12 @@
integer, parameter :: NGLLSQUARE_NDIM = NGLLSQUARE * NDIM
integer, parameter :: NGLLCUBE_NDIM = NGLLCUBE * NDIM
-! number of nodes for 2D and 3D shape functions for hexahedra
-! we use 8-node mesh bricks, which are more stable than 27-node elements
- integer, parameter :: NGNOD = 8, NGNOD2D = 4
+! this to only take the corners (i.e. extract a QUAD4 or HEX8 in some routines even if we use higher-order elements)
integer, parameter :: NGNOD_EIGHT_CORNERS = 8
+ integer, parameter :: NGNOD2D_FOUR_CORNERS = 4
! number of points in each AVS or OpenDX quadrangular cell for movies
- integer, parameter :: NGNOD2D_AVS_DX = 4
+ integer, parameter :: NGNOD2D_FOUR_CORNERS_AVS_DX = NGNOD2D_FOUR_CORNERS
! a few useful constants
double precision, parameter :: ZERO = 0.d0,ONE = 1.d0,TWO = 2.d0 !,HALF = 0.5d0
Modified: seismo/3D/SPECFEM3D/trunk/src/shared/create_movie_shakemap_AVS_DX_GMT.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/shared/create_movie_shakemap_AVS_DX_GMT.f90 2012-08-20 12:15:45 UTC (rev 20598)
+++ seismo/3D/SPECFEM3D/trunk/src/shared/create_movie_shakemap_AVS_DX_GMT.f90 2012-08-20 12:57:46 UTC (rev 20599)
@@ -113,7 +113,7 @@
!!!! NL NL
! order of points representing the 2D square element
- integer,dimension(NGNOD2D_AVS_DX),parameter :: iorder = (/1,3,2,4/)
+ integer,dimension(NGNOD2D_FOUR_CORNERS_AVS_DX),parameter :: iorder = (/1,3,2,4/)
! ************** PROGRAM STARTS HERE **************
@@ -145,7 +145,7 @@
if(USE_HIGHRES_FOR_MOVIES) then
ilocnum = NSPEC_SURFACE_EXT_MESH*NGLLSQUARE
else
- ilocnum = NSPEC_SURFACE_EXT_MESH*NGNOD2D_AVS_DX
+ ilocnum = NSPEC_SURFACE_EXT_MESH*NGNOD2D_FOUR_CORNERS_AVS_DX
endif
print*,' moviedata element surfaces: ',NSPEC_SURFACE_EXT_MESH
print*,' moviedata total elements all: ',ilocnum
@@ -250,7 +250,7 @@
endif
! maximum theoretical number of points at the surface
- npointot = NGNOD2D_AVS_DX * nspectot_AVS_max
+ npointot = NGNOD2D_FOUR_CORNERS_AVS_DX * nspectot_AVS_max
! allocate arrays for sorting routine
allocate(iglob(npointot),loc(npointot), &
@@ -407,8 +407,8 @@
do j = 1,NGLLY-1
do i = 1,NGLLX-1
- ieoff = NGNOD2D_AVS_DX*(ielm+(i-1)+(j-1)*(NGLLX-1))
- do ilocnum = 1,NGNOD2D_AVS_DX
+ ieoff = NGNOD2D_FOUR_CORNERS_AVS_DX*(ielm+(i-1)+(j-1)*(NGLLX-1))
+ do ilocnum = 1,NGNOD2D_FOUR_CORNERS_AVS_DX
if(ilocnum == 1) then
xp(ieoff+ilocnum) = dble(x(i,j))
@@ -446,10 +446,10 @@
else
! low-resolution (only spectral element corners)
ispec = ispec + 1
- ieoff = NGNOD2D_AVS_DX*(ispec-1)
+ ieoff = NGNOD2D_FOUR_CORNERS_AVS_DX*(ispec-1)
! four points for each element
- do i = 1,NGNOD2D_AVS_DX
+ do i = 1,NGNOD2D_FOUR_CORNERS_AVS_DX
! accounts for different ordering of square points
ilocnum = iorder(i)
@@ -657,9 +657,9 @@
! output list of points
mask_point = .false.
do ispec=1,nspectot_AVS_max
- ieoff = NGNOD2D_AVS_DX*(ispec-1)
+ ieoff = NGNOD2D_FOUR_CORNERS_AVS_DX*(ispec-1)
! four points for each element
- do ilocnum = 1,NGNOD2D_AVS_DX
+ do ilocnum = 1,NGNOD2D_FOUR_CORNERS_AVS_DX
ibool_number = iglob(ilocnum+ieoff)
if(.not. mask_point(ibool_number)) then
call utm_geo(long,lat,xp_save(ilocnum+ieoff),yp_save(ilocnum+ieoff), &
@@ -676,9 +676,9 @@
mask_point = .false.
ipoin = 0
do ispec=1,nspectot_AVS_max
- ieoff = NGNOD2D_AVS_DX*(ispec-1)
+ ieoff = NGNOD2D_FOUR_CORNERS_AVS_DX*(ispec-1)
! four points for each element
- do ilocnum = 1,NGNOD2D_AVS_DX
+ do ilocnum = 1,NGNOD2D_FOUR_CORNERS_AVS_DX
ibool_number = iglob(ilocnum+ieoff)
if(.not. mask_point(ibool_number)) then
ipoin = ipoin + 1
@@ -699,7 +699,7 @@
! output list of elements
do ispec=1,nspectot_AVS_max
- ieoff = NGNOD2D_AVS_DX*(ispec-1)
+ ieoff = NGNOD2D_FOUR_CORNERS_AVS_DX*(ispec-1)
! four points for each element
ibool_number1 = iglob(ieoff + 1)
ibool_number2 = iglob(ieoff + 2)
@@ -728,9 +728,9 @@
! output data values
mask_point = .false.
do ispec=1,nspectot_AVS_max
- ieoff = NGNOD2D_AVS_DX*(ispec-1)
+ ieoff = NGNOD2D_FOUR_CORNERS_AVS_DX*(ispec-1)
! four points for each element
- do ilocnum = 1,NGNOD2D_AVS_DX
+ do ilocnum = 1,NGNOD2D_FOUR_CORNERS_AVS_DX
ibool_number = iglob(ilocnum+ieoff)
if(.not. mask_point(ibool_number)) then
if(USE_OPENDX) then
@@ -851,8 +851,8 @@
! establish initial pointers
do ispec=1,nspec
- ieoff=NGNOD2D_AVS_DX*(ispec-1)
- do ilocnum=1,NGNOD2D_AVS_DX
+ ieoff=NGNOD2D_FOUR_CORNERS_AVS_DX*(ispec-1)
+ do ilocnum=1,NGNOD2D_FOUR_CORNERS_AVS_DX
loc(ilocnum+ieoff)=ilocnum+ieoff
enddo
enddo
Modified: seismo/3D/SPECFEM3D/trunk/src/specfem3D/setup_movie_meshes.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/specfem3D/setup_movie_meshes.f90 2012-08-20 12:15:45 UTC (rev 20598)
+++ seismo/3D/SPECFEM3D/trunk/src/specfem3D/setup_movie_meshes.f90 2012-08-20 12:57:46 UTC (rev 20599)
@@ -56,13 +56,13 @@
store_val_uz_external_mesh(NGLLX*NGLLY*1),stat=ier)
if( ier /= 0 ) stop 'error allocating dummy arrays for highres movie'
else
- allocate(faces_surface_ext_mesh(NGNOD2D,1), &
- store_val_x_external_mesh(NGNOD2D*1), &
- store_val_y_external_mesh(NGNOD2D*1), &
- store_val_z_external_mesh(NGNOD2D*1), &
- store_val_ux_external_mesh(NGNOD2D*1), &
- store_val_uy_external_mesh(NGNOD2D*1), &
- store_val_uz_external_mesh(NGNOD2D*1),stat=ier)
+ allocate(faces_surface_ext_mesh(NGNOD2D_FOUR_CORNERS,1), &
+ store_val_x_external_mesh(NGNOD2D_FOUR_CORNERS*1), &
+ store_val_y_external_mesh(NGNOD2D_FOUR_CORNERS*1), &
+ store_val_z_external_mesh(NGNOD2D_FOUR_CORNERS*1), &
+ store_val_ux_external_mesh(NGNOD2D_FOUR_CORNERS*1), &
+ store_val_uy_external_mesh(NGNOD2D_FOUR_CORNERS*1), &
+ store_val_uz_external_mesh(NGNOD2D_FOUR_CORNERS*1),stat=ier)
if( ier /= 0 ) stop 'error allocating dummy arrays for lowres movie'
endif
else
@@ -76,13 +76,13 @@
store_val_uz_external_mesh(NGLLX*NGLLY*nfaces_surface_ext_mesh),stat=ier)
if( ier /= 0 ) stop 'error allocating arrays for highres movie'
else
- allocate(faces_surface_ext_mesh(NGNOD2D,nfaces_surface_ext_mesh), &
- store_val_x_external_mesh(NGNOD2D*nfaces_surface_ext_mesh), &
- store_val_y_external_mesh(NGNOD2D*nfaces_surface_ext_mesh), &
- store_val_z_external_mesh(NGNOD2D*nfaces_surface_ext_mesh), &
- store_val_ux_external_mesh(NGNOD2D*nfaces_surface_ext_mesh), &
- store_val_uy_external_mesh(NGNOD2D*nfaces_surface_ext_mesh), &
- store_val_uz_external_mesh(NGNOD2D*nfaces_surface_ext_mesh),stat=ier)
+ allocate(faces_surface_ext_mesh(NGNOD2D_FOUR_CORNERS,nfaces_surface_ext_mesh), &
+ store_val_x_external_mesh(NGNOD2D_FOUR_CORNERS*nfaces_surface_ext_mesh), &
+ store_val_y_external_mesh(NGNOD2D_FOUR_CORNERS*nfaces_surface_ext_mesh), &
+ store_val_z_external_mesh(NGNOD2D_FOUR_CORNERS*nfaces_surface_ext_mesh), &
+ store_val_ux_external_mesh(NGNOD2D_FOUR_CORNERS*nfaces_surface_ext_mesh), &
+ store_val_uy_external_mesh(NGNOD2D_FOUR_CORNERS*nfaces_surface_ext_mesh), &
+ store_val_uz_external_mesh(NGNOD2D_FOUR_CORNERS*nfaces_surface_ext_mesh),stat=ier)
if( ier /= 0 ) stop 'error allocating arrays for lowres movie'
endif
endif
@@ -104,12 +104,12 @@
store_val_uz_all_external_mesh(NGLLX*NGLLY*nfaces_surface_glob_ext_mesh),stat=ier)
if( ier /= 0 ) stop 'error allocating arrays for highres movie'
else
- allocate(store_val_x_all_external_mesh(NGNOD2D*nfaces_surface_glob_ext_mesh), &
- store_val_y_all_external_mesh(NGNOD2D*nfaces_surface_glob_ext_mesh), &
- store_val_z_all_external_mesh(NGNOD2D*nfaces_surface_glob_ext_mesh), &
- store_val_ux_all_external_mesh(NGNOD2D*nfaces_surface_glob_ext_mesh), &
- store_val_uy_all_external_mesh(NGNOD2D*nfaces_surface_glob_ext_mesh), &
- store_val_uz_all_external_mesh(NGNOD2D*nfaces_surface_glob_ext_mesh),stat=ier)
+ allocate(store_val_x_all_external_mesh(NGNOD2D_FOUR_CORNERS*nfaces_surface_glob_ext_mesh), &
+ store_val_y_all_external_mesh(NGNOD2D_FOUR_CORNERS*nfaces_surface_glob_ext_mesh), &
+ store_val_z_all_external_mesh(NGNOD2D_FOUR_CORNERS*nfaces_surface_glob_ext_mesh), &
+ store_val_ux_all_external_mesh(NGNOD2D_FOUR_CORNERS*nfaces_surface_glob_ext_mesh), &
+ store_val_uy_all_external_mesh(NGNOD2D_FOUR_CORNERS*nfaces_surface_glob_ext_mesh), &
+ store_val_uz_all_external_mesh(NGNOD2D_FOUR_CORNERS*nfaces_surface_glob_ext_mesh),stat=ier)
if( ier /= 0 ) stop 'error allocating arrays for lowres movie'
endif
endif
@@ -123,7 +123,7 @@
if (USE_HIGHRES_FOR_MOVIES) then
faces_surface_offset_ext_mesh(:) = faces_surface_offset_ext_mesh(:)*NGLLX*NGLLY
else
- faces_surface_offset_ext_mesh(:) = faces_surface_offset_ext_mesh(:)*NGNOD2D
+ faces_surface_offset_ext_mesh(:) = faces_surface_offset_ext_mesh(:)*NGNOD2D_FOUR_CORNERS
endif
! stores global indices of GLL points on the surface to array faces_surface_ext_mesh
@@ -304,9 +304,9 @@
nfaces_surface_glob_em_points = nfaces_surface_glob_ext_mesh*NGLLX*NGLLY
else
! low-res movies only output at element corners
- nfaces_perproc_surface_ext_mesh(:) = nfaces_perproc_surface_ext_mesh(:)*NGNOD2D
- nfaces_surface_ext_mesh_points = nfaces_surface_ext_mesh*NGNOD2D
- nfaces_surface_glob_em_points = nfaces_surface_glob_ext_mesh*NGNOD2D
+ nfaces_perproc_surface_ext_mesh(:) = nfaces_perproc_surface_ext_mesh(:)*NGNOD2D_FOUR_CORNERS
+ nfaces_surface_ext_mesh_points = nfaces_surface_ext_mesh*NGNOD2D_FOUR_CORNERS
+ nfaces_surface_glob_em_points = nfaces_surface_glob_ext_mesh*NGNOD2D_FOUR_CORNERS
endif
end subroutine setup_movie_meshes
Modified: seismo/3D/SPECFEM3D/trunk/src/specfem3D/specfem3D_par.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/specfem3D/specfem3D_par.f90 2012-08-20 12:15:45 UTC (rev 20598)
+++ seismo/3D/SPECFEM3D/trunk/src/specfem3D/specfem3D_par.f90 2012-08-20 12:57:46 UTC (rev 20599)
@@ -518,7 +518,7 @@
! parameter module for movies/shakemovies
- use constants,only: CUSTOM_REAL,NGLLX,NGLLY,NGLLZ,NGNOD2D
+ use constants,only: CUSTOM_REAL,NGLLX,NGLLY,NGLLZ,NGNOD2D_FOUR_CORNERS
implicit none
@@ -557,7 +557,7 @@
integer :: nfaces_surface_ext_mesh,nfaces_surface_ext_mesh_points
integer :: nfaces_surface_glob_ext_mesh,nfaces_surface_glob_em_points
! face corner indices
- integer :: iorderi(NGNOD2D),iorderj(NGNOD2D)
+ integer :: iorderi(NGNOD2D_FOUR_CORNERS),iorderj(NGNOD2D_FOUR_CORNERS)
! movie parameters
double precision :: HDUR_MOVIE
Modified: seismo/3D/SPECFEM3D/trunk/src/specfem3D/write_movie_output.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/specfem3D/write_movie_output.f90 2012-08-20 12:15:45 UTC (rev 20598)
+++ seismo/3D/SPECFEM3D/trunk/src/specfem3D/write_movie_output.f90 2012-08-20 12:57:46 UTC (rev 20599)
@@ -131,9 +131,9 @@
do ipoin = 1, 4
iglob = faces_surface_ext_mesh(ipoin,ispec2D)
! x,y,z coordinates
- store_val_x_external_mesh(NGNOD2D*(ispec2D-1)+ipoin) = xstore(iglob)
- store_val_y_external_mesh(NGNOD2D*(ispec2D-1)+ipoin) = ystore(iglob)
- store_val_z_external_mesh(NGNOD2D*(ispec2D-1)+ipoin) = zstore(iglob)
+ store_val_x_external_mesh(NGNOD2D_FOUR_CORNERS*(ispec2D-1)+ipoin) = xstore(iglob)
+ store_val_y_external_mesh(NGNOD2D_FOUR_CORNERS*(ispec2D-1)+ipoin) = ystore(iglob)
+ store_val_z_external_mesh(NGNOD2D_FOUR_CORNERS*(ispec2D-1)+ipoin) = zstore(iglob)
enddo
endif
enddo
@@ -195,22 +195,22 @@
enddo
else
- ! low-resolution: only corner points outputted
+ ! low-resolution: only corner points are output
do ipoin = 1, 4
iglob = faces_surface_ext_mesh(ipoin,ispec2D)
! saves norm of displacement,velocity and acceleration vector
if( ispec_is_elastic(ispec) ) then
! norm of displacement
- store_val_ux_external_mesh(NGNOD2D*(ispec2D-1)+ipoin) = &
- max(store_val_ux_external_mesh(NGNOD2D*(ispec2D-1)+ipoin), &
+ store_val_ux_external_mesh(NGNOD2D_FOUR_CORNERS*(ispec2D-1)+ipoin) = &
+ max(store_val_ux_external_mesh(NGNOD2D_FOUR_CORNERS*(ispec2D-1)+ipoin), &
sqrt(displ(1,iglob)**2 + displ(2,iglob)**2 + displ(3,iglob)**2))
! norm of velocity
- store_val_uy_external_mesh(NGNOD2D*(ispec2D-1)+ipoin) = &
- max(store_val_uy_external_mesh(NGNOD2D*(ispec2D-1)+ipoin), &
+ store_val_uy_external_mesh(NGNOD2D_FOUR_CORNERS*(ispec2D-1)+ipoin) = &
+ max(store_val_uy_external_mesh(NGNOD2D_FOUR_CORNERS*(ispec2D-1)+ipoin), &
sqrt(veloc(1,iglob)**2 + veloc(2,iglob)**2 + veloc(3,iglob)**2))
! norm of acceleration
- store_val_uz_external_mesh(NGNOD2D*(ispec2D-1)+ipoin) = &
- max(store_val_uz_external_mesh(NGNOD2D*(ispec2D-1)+ipoin), &
+ store_val_uz_external_mesh(NGNOD2D_FOUR_CORNERS*(ispec2D-1)+ipoin) = &
+ max(store_val_uz_external_mesh(NGNOD2D_FOUR_CORNERS*(ispec2D-1)+ipoin), &
sqrt(accel(1,iglob)**2 + accel(2,iglob)**2 + accel(3,iglob)**2))
endif
@@ -219,7 +219,7 @@
! sets velocity vector with maximum norm of wavefield values
call wmo_get_max_vector(ispec,ispec2D,iglob,ipoin, &
displ_element,veloc_element,accel_element, &
- NGNOD2D)
+ NGNOD2D_FOUR_CORNERS)
endif
enddo
endif
@@ -318,20 +318,20 @@
store_val_ux_external_mesh(narraydim*(ispec2D-1)+ipoin) = &
max(store_val_ux_external_mesh(narraydim*(ispec2D-1)+ipoin), &
sqrt(displ_element(1,i,j,k)**2 &
- + displ_element(2,i,j,k)**2 &
- + displ_element(3,i,j,k)**2))
+ + displ_element(2,i,j,k)**2 &
+ + displ_element(3,i,j,k)**2))
! norm of velocity
store_val_uy_external_mesh(narraydim*(ispec2D-1)+ipoin) = &
max(store_val_uy_external_mesh(narraydim*(ispec2D-1)+ipoin), &
sqrt(veloc_element(1,i,j,k)**2 &
- + veloc_element(2,i,j,k)**2 &
- + veloc_element(3,i,j,k)**2))
+ + veloc_element(2,i,j,k)**2 &
+ + veloc_element(3,i,j,k)**2))
! norm of acceleration
store_val_uz_external_mesh(narraydim*(ispec2D-1)+ipoin) = &
max(store_val_uz_external_mesh(narraydim*(ispec2D-1)+ipoin), &
sqrt(accel_element(1,i,j,k)**2 &
- + accel_element(2,i,j,k)**2 &
- + accel_element(3,i,j,k)**2))
+ + accel_element(2,i,j,k)**2 &
+ + accel_element(3,i,j,k)**2))
! not really needed, but could be used to check...
is_done = .true.
return
@@ -377,9 +377,9 @@
do ipoin = 1, 4
iglob = faces_surface_ext_mesh(ipoin,ispec2D)
! x,y,z coordinates
- store_val_x_external_mesh(NGNOD2D*(ispec2D-1)+ipoin) = xstore(iglob)
- store_val_y_external_mesh(NGNOD2D*(ispec2D-1)+ipoin) = ystore(iglob)
- store_val_z_external_mesh(NGNOD2D*(ispec2D-1)+ipoin) = zstore(iglob)
+ store_val_x_external_mesh(NGNOD2D_FOUR_CORNERS*(ispec2D-1)+ipoin) = xstore(iglob)
+ store_val_y_external_mesh(NGNOD2D_FOUR_CORNERS*(ispec2D-1)+ipoin) = ystore(iglob)
+ store_val_z_external_mesh(NGNOD2D_FOUR_CORNERS*(ispec2D-1)+ipoin) = zstore(iglob)
enddo
endif
enddo
@@ -423,7 +423,7 @@
! puts displ/velocity values into storage array
call wmo_get_vel_vector(ispec,ispec2D,ipoin,iglob, &
val_element, &
- NGNOD2D)
+ NGNOD2D_FOUR_CORNERS)
enddo
endif
enddo
@@ -603,7 +603,7 @@
jmax = maxval( free_surface_ijk(2,:,iface) )
kmin = minval( free_surface_ijk(3,:,iface) )
kmax = maxval( free_surface_ijk(3,:,iface) )
- do iloc = 1, NGNOD2D
+ do iloc = 1, NGNOD2D_FOUR_CORNERS
ipoin = ipoin + 1
! corner points
if( imin == imax ) then
@@ -669,7 +669,7 @@
jmax = maxval( free_surface_ijk(2,:,iface) )
kmin = minval( free_surface_ijk(3,:,iface) )
kmax = maxval( free_surface_ijk(3,:,iface) )
- do iloc = 1, NGNOD2D
+ do iloc = 1, NGNOD2D_FOUR_CORNERS
ipoin = ipoin + 1
! corner points
if( imin == imax ) then
@@ -845,7 +845,7 @@
jmax = maxval( free_surface_ijk(2,:,iface) )
kmin = minval( free_surface_ijk(3,:,iface) )
kmax = maxval( free_surface_ijk(3,:,iface) )
- do iloc = 1, NGNOD2D
+ do iloc = 1, NGNOD2D_FOUR_CORNERS
ipoin = ipoin + 1
! corner points
if( imin == imax ) then
More information about the CIG-COMMITS
mailing list