[cig-commits] r20837 - in seismo/3D/SPECFEM3D/trunk/src: decompose_mesh_SCOTCH generate_databases shared specfem3D
dkomati1 at geodynamics.org
dkomati1 at geodynamics.org
Sat Oct 13 15:18:10 PDT 2012
Author: dkomati1
Date: 2012-10-13 15:18:10 -0700 (Sat, 13 Oct 2012)
New Revision: 20837
Modified:
seismo/3D/SPECFEM3D/trunk/src/decompose_mesh_SCOTCH/decompose_mesh_SCOTCH.f90
seismo/3D/SPECFEM3D/trunk/src/decompose_mesh_SCOTCH/part_decompose_mesh_SCOTCH.f90
seismo/3D/SPECFEM3D/trunk/src/generate_databases/create_regions_mesh.f90
seismo/3D/SPECFEM3D/trunk/src/generate_databases/generate_databases.f90
seismo/3D/SPECFEM3D/trunk/src/generate_databases/get_absorbing_boundary.f90
seismo/3D/SPECFEM3D/trunk/src/generate_databases/get_coupling_surfaces.f90
seismo/3D/SPECFEM3D/trunk/src/generate_databases/read_partition_files.f90
seismo/3D/SPECFEM3D/trunk/src/shared/get_element_face.f90
seismo/3D/SPECFEM3D/trunk/src/shared/get_force.f90
seismo/3D/SPECFEM3D/trunk/src/shared/get_jacobian_boundaries.f90
seismo/3D/SPECFEM3D/trunk/src/shared/get_shape2D.f90
seismo/3D/SPECFEM3D/trunk/src/shared/get_shape3D.f90
seismo/3D/SPECFEM3D/trunk/src/shared/prepare_assemble_MPI.f90
seismo/3D/SPECFEM3D/trunk/src/shared/read_parameter_file.f90
seismo/3D/SPECFEM3D/trunk/src/shared/save_header_file.f90
seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_add_sources_elastic.f90
seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_add_sources_poroelastic.f90
seismo/3D/SPECFEM3D/trunk/src/specfem3D/initialize_simulation.f90
seismo/3D/SPECFEM3D/trunk/src/specfem3D/locate_source.f90
seismo/3D/SPECFEM3D/trunk/src/specfem3D/model_update.f90
seismo/3D/SPECFEM3D/trunk/src/specfem3D/noise_tomography.f90
seismo/3D/SPECFEM3D/trunk/src/specfem3D/prepare_timerun.F90
seismo/3D/SPECFEM3D/trunk/src/specfem3D/save_external_bin_m_up.f90
seismo/3D/SPECFEM3D/trunk/src/specfem3D/setup_sources_receivers.f90
Log:
done adding support for 27-node elements (not tested yet);
also removed useless white spaces
Modified: seismo/3D/SPECFEM3D/trunk/src/decompose_mesh_SCOTCH/decompose_mesh_SCOTCH.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/decompose_mesh_SCOTCH/decompose_mesh_SCOTCH.f90 2012-10-13 18:42:30 UTC (rev 20836)
+++ seismo/3D/SPECFEM3D/trunk/src/decompose_mesh_SCOTCH/decompose_mesh_SCOTCH.f90 2012-10-13 22:18:10 UTC (rev 20837)
@@ -24,7 +24,6 @@
!
!=====================================================================
-
module decompose_mesh_SCOTCH
use part_decompose_mesh_SCOTCH
@@ -62,12 +61,12 @@
integer :: ninterfaces
integer :: my_ninterface
- integer :: nsize ! Max number of elements that contain the same node.
+ integer :: nsize ! max number of elements that contain the same node
integer :: nb_edges
integer :: ispec, inode
- integer :: max_neighbour ! Real maximum number of neighbours per element
- integer :: sup_neighbour ! Majoration of the maximum number of neighbours per element
+ integer :: max_neighbour ! real maximum number of neighbours per element
+ integer :: sup_neighbour ! majoration (overestimate) of the maximum number of neighbours per element
integer :: ipart, nnodes_loc, nspec_local
integer :: num_elmnt, num_node, num_mat
@@ -116,10 +115,13 @@
! reads in mesh files
!----------------------------------------------------------------------------------------------
subroutine read_mesh_files
+
implicit none
+
character(len=256) :: line
logical :: use_poroelastic_file
integer(long) :: nspec_long
+ integer :: inode
! reads node coordinates
open(unit=98, file=localpath_name(1:len_trim(localpath_name))//'/nodes_coords_file',&
@@ -134,7 +136,7 @@
do inode = 1, nnodes
! format: #id_node #x_coordinate #y_coordinate #z_coordinate
read(98,*) num_node, nodes_coords(1,num_node), nodes_coords(2,num_node), nodes_coords(3,num_node)
- end do
+ enddo
close(98)
print*, 'total number of nodes: '
print*, ' nnodes = ', nnodes
@@ -148,7 +150,7 @@
read(98,*) nspec_long
! debug check size limit
- if( nspec_long > 2147483647 ) then
+ if( nspec_long > 2147483646 ) then
print *,'size exceeds integer 4-byte limit: ',nspec_long
print*,'bit size fortran: ',bit_size(nspec)
stop 'error number of elements too large'
@@ -161,6 +163,7 @@
if( ier /= 0 ) stop 'error allocating array elmnts'
do ispec = 1, nspec
! format: # element_id #id_node1 ... #id_node8
+ ! or # element_id #id_node1 ... #id_node27
! note: be aware that here we can have different node ordering for a cube element;
! the ordering from Cubit files might not be consistent for multiple volumes, or uneven, unstructured grids
@@ -171,12 +174,11 @@
! then top (positive z-direction) of element
! point 5 = (0,0,1), point 6 = (0,1,1), point 7 = (1,1,1), point 8 = (1,0,1)
- read(98,*) num_elmnt, elmnts(1,num_elmnt), elmnts(2,num_elmnt),elmnts(3,num_elmnt), elmnts(4,num_elmnt), &
- elmnts(5,num_elmnt), elmnts(6,num_elmnt), elmnts(7,num_elmnt), elmnts(8,num_elmnt)
+ read(98,*) num_elmnt,(elmnts(inode,num_elmnt), inode=1,NGNOD)
if((num_elmnt > nspec) .or. (num_elmnt < 1) ) stop "ERROR : Invalid mesh file."
- end do
+ enddo
close(98)
print*, 'total number of spectral elements:'
print*, ' nspec = ', nspec
@@ -193,7 +195,7 @@
! note: be aware that elements may not be sorted in materials_file
read(98,*) num_mat,mat(1,num_mat)
if((num_mat > nspec) .or. (num_mat < 1) ) stop "ERROR : Invalid mat file."
- end do
+ enddo
close(98)
! TODO:
@@ -235,9 +237,9 @@
else
! negative materials_id: undefined material properties yet
count_undef_mat = count_undef_mat + 1
- end if
+ endif
read(98,*,iostat=ier) idummy,num_mat
- end do
+ enddo
close(98)
print*, ' defined = ',count_def_mat, 'undefined = ',count_undef_mat
! check with material flags
@@ -341,7 +343,7 @@
endif !if(idomain_id == 1 .or. idomain_id == 2)
- end do
+ enddo
! reads in undefined material properties
rewind(98,iostat=ier) ! back to the beginning of the file
@@ -410,7 +412,7 @@
stop "ERROR: invalid flag_up in interface definition in nummaterial_velocity_file"
endif
endif
- end do
+ enddo
if( use_poroelastic_file ) close(97)
close(98)
@@ -452,7 +454,7 @@
endif
allocate(ibelm_xmin(nspec2D_xmin),stat=ier)
if( ier /= 0 ) stop 'error allocating array ibelm_xmin'
- allocate(nodes_ibelm_xmin(4,nspec2D_xmin),stat=ier)
+ allocate(nodes_ibelm_xmin(NGNOD2D,nspec2D_xmin),stat=ier)
if( ier /= 0 ) stop 'error allocating array nodes_ibelm_xmin'
do ispec2D = 1,nspec2D_xmin
! format: #id_(element containing the face) #id_node1_face .. #id_node4_face
@@ -462,10 +464,8 @@
!
! doesn't necessarily have to start on top-rear, then bottom-rear, bottom-front, and finally top-front i.e.:
! point 1 = (0,1,1), point 2 = (0,1,0), point 3 = (0,0,0), point 4 = (0,0,1)
- read(98,*) ibelm_xmin(ispec2D), nodes_ibelm_xmin(1,ispec2D), nodes_ibelm_xmin(2,ispec2D), &
- nodes_ibelm_xmin(3,ispec2D), nodes_ibelm_xmin(4,ispec2D) !! DK DK see if we need NGNOD here
-
- end do
+ read(98,*) ibelm_xmin(ispec2D), (nodes_ibelm_xmin(inode,ispec2D), inode=1,NGNOD2D)
+ enddo
close(98)
print*, 'absorbing boundaries:'
print*, ' nspec2D_xmin = ', nspec2D_xmin
@@ -480,13 +480,12 @@
endif
allocate(ibelm_xmax(nspec2D_xmax),stat=ier)
if( ier /= 0 ) stop 'error allocating array ibelm_xmax'
- allocate(nodes_ibelm_xmax(4,nspec2D_xmax),stat=ier)
+ allocate(nodes_ibelm_xmax(NGNOD2D,nspec2D_xmax),stat=ier)
if( ier /= 0 ) stop 'error allocating array nodes_ibelm_xmax'
do ispec2D = 1,nspec2D_xmax
! format: #id_(element containing the face) #id_node1_face .. #id_node4_face
- read(98,*) ibelm_xmax(ispec2D), nodes_ibelm_xmax(1,ispec2D), nodes_ibelm_xmax(2,ispec2D), &
- nodes_ibelm_xmax(3,ispec2D), nodes_ibelm_xmax(4,ispec2D)
- end do
+ read(98,*) ibelm_xmax(ispec2D), (nodes_ibelm_xmax(inode,ispec2D), inode=1,NGNOD2D)
+ enddo
close(98)
print*, ' nspec2D_xmax = ', nspec2D_xmax
@@ -500,13 +499,12 @@
endif
allocate(ibelm_ymin(nspec2D_ymin),stat=ier)
if( ier /= 0 ) stop 'error allocating array ibelm_ymin'
- allocate(nodes_ibelm_ymin(4,nspec2D_ymin),stat=ier)
+ allocate(nodes_ibelm_ymin(NGNOD2D,nspec2D_ymin),stat=ier)
if( ier /= 0 ) stop 'error allocating array nodes_ibelm_ymin'
do ispec2D = 1,nspec2D_ymin
! format: #id_(element containing the face) #id_node1_face .. #id_node4_face
- read(98,*) ibelm_ymin(ispec2D), nodes_ibelm_ymin(1,ispec2D), nodes_ibelm_ymin(2,ispec2D), &
- nodes_ibelm_ymin(3,ispec2D), nodes_ibelm_ymin(4,ispec2D)
- end do
+ read(98,*) ibelm_ymin(ispec2D), (nodes_ibelm_ymin(inode,ispec2D), inode=1,NGNOD2D)
+ enddo
close(98)
print*, ' nspec2D_ymin = ', nspec2D_ymin
@@ -520,13 +518,12 @@
endif
allocate(ibelm_ymax(nspec2D_ymax),stat=ier)
if( ier /= 0 ) stop 'error allocating array ibelm_ymax'
- allocate(nodes_ibelm_ymax(4,nspec2D_ymax),stat=ier)
+ allocate(nodes_ibelm_ymax(NGNOD2D,nspec2D_ymax),stat=ier)
if( ier /= 0 ) stop 'error allocating array nodes_ibelm_ymax'
do ispec2D = 1,nspec2D_ymax
! format: #id_(element containing the face) #id_node1_face .. #id_node4_face
- read(98,*) ibelm_ymax(ispec2D), nodes_ibelm_ymax(1,ispec2D), nodes_ibelm_ymax(2,ispec2D), &
- nodes_ibelm_ymax(3,ispec2D), nodes_ibelm_ymax(4,ispec2D)
- end do
+ read(98,*) ibelm_ymax(ispec2D), (nodes_ibelm_ymax(inode,ispec2D), inode=1,NGNOD2D)
+ enddo
close(98)
print*, ' nspec2D_ymax = ', nspec2D_ymax
@@ -540,13 +537,12 @@
endif
allocate(ibelm_bottom(nspec2D_bottom),stat=ier)
if( ier /= 0 ) stop 'error allocating array ibelm_bottom'
- allocate(nodes_ibelm_bottom(4,nspec2D_bottom),stat=ier)
+ allocate(nodes_ibelm_bottom(NGNOD2D,nspec2D_bottom),stat=ier)
if( ier /= 0 ) stop 'error allocating array nodes_ibelm_bottom'
do ispec2D = 1,nspec2D_bottom
! format: #id_(element containing the face) #id_node1_face .. #id_node4_face
- read(98,*) ibelm_bottom(ispec2D), nodes_ibelm_bottom(1,ispec2D), nodes_ibelm_bottom(2,ispec2D), &
- nodes_ibelm_bottom(3,ispec2D), nodes_ibelm_bottom(4,ispec2D)
- end do
+ read(98,*) ibelm_bottom(ispec2D), (nodes_ibelm_bottom(inode,ispec2D), inode=1,NGNOD2D)
+ enddo
close(98)
print*, ' nspec2D_bottom = ', nspec2D_bottom
@@ -560,13 +556,12 @@
endif
allocate(ibelm_top(nspec2D_top),stat=ier)
if( ier /= 0 ) stop 'error allocating array ibelm_top'
- allocate(nodes_ibelm_top(4,nspec2D_top),stat=ier)
+ allocate(nodes_ibelm_top(NGNOD2D,nspec2D_top),stat=ier)
if( ier /= 0 ) stop 'error allocating array nodes_ibelm_top'
do ispec2D = 1,nspec2D_top
! format: #id_(element containing the face) #id_node1_face .. #id_node4_face
- read(98,*) ibelm_top(ispec2D), nodes_ibelm_top(1,ispec2D), nodes_ibelm_top(2,ispec2D), &
- nodes_ibelm_top(3,ispec2D), nodes_ibelm_top(4,ispec2D)
- end do
+ read(98,*) ibelm_top(ispec2D), (nodes_ibelm_top(inode,ispec2D), inode=1,NGNOD2D)
+ enddo
close(98)
print*, ' nspec2D_top = ', nspec2D_top
@@ -580,13 +575,12 @@
endif
allocate(ibelm_moho(nspec2D_moho),stat=ier)
if( ier /= 0 ) stop 'error allocating array ibelm_moho'
- allocate(nodes_ibelm_moho(4,nspec2D_moho),stat=ier)
+ allocate(nodes_ibelm_moho(NGNOD2D,nspec2D_moho),stat=ier)
if( ier /= 0 ) stop 'error allocating array nodes_ibelm_moho'
do ispec2D = 1,nspec2D_moho
! format: #id_(element containing the face) #id_node1_face .. #id_node4_face
- read(98,*) ibelm_moho(ispec2D), nodes_ibelm_moho(1,ispec2D), nodes_ibelm_moho(2,ispec2D), &
- nodes_ibelm_moho(3,ispec2D), nodes_ibelm_moho(4,ispec2D)
- end do
+ read(98,*) ibelm_moho(ispec2D), (nodes_ibelm_moho(inode,ispec2D), inode=1,NGNOD2D)
+ enddo
close(98)
if( nspec2D_moho > 0 ) print*, ' nspec2D_moho = ', nspec2D_moho
@@ -605,30 +599,38 @@
mask_nodes_elmnts(:) = .false.
used_nodes_elmnts(:) = 0
do ispec = 1, nspec
- do inode = 1, NGNOD
+ do inode = 1, NGNOD_EIGHT_CORNERS
mask_nodes_elmnts(elmnts(inode,ispec)) = .true.
used_nodes_elmnts(elmnts(inode,ispec)) = used_nodes_elmnts(elmnts(inode,ispec)) + 1
enddo
enddo
- print *, 'nodes valence: '
- print *, ' min = ',minval(used_nodes_elmnts(:)),'max = ', maxval(used_nodes_elmnts(:))
+ print *, 'node valence:'
+ print *, ' min = ',minval(used_nodes_elmnts(:)),' max = ', maxval(used_nodes_elmnts(:))
do inode = 1, nnodes
if (.not. mask_nodes_elmnts(inode)) then
- stop 'ERROR : nodes not used.'
+ stop 'ERROR: found some unused nodes (weird, but not necessarily fatal; your mesher may have created extra nodes).'
endif
enddo
+
+! max number of elements that contain the same node
nsize = maxval(used_nodes_elmnts(:))
+! majoration (overestimate) of the maximum number of neighbours per element
+!! DK DK nfaces is a constant equal to 6 (number of faces of a cube).
+!! DK DK I have no idea how this formula works; it was designed by Nicolas Le Goff
+ sup_neighbour = NGNOD_EIGHT_CORNERS * nsize - (NGNOD_EIGHT_CORNERS + (NGNOD_EIGHT_CORNERS/2 - 1)*nfaces)
+
! debug check size limit
- if( NGNOD * nsize - (NGNOD + (NGNOD/2 - 1)*nfaces) > 2147483647 ) then
- print *,'size exceeds integer 4-byte limit: ',sup_neighbour,NGNOD,nsize,nfaces
- print*,'bit size fortran: ',bit_size(sup_neighbour)
+!! DK DK this check will likely fail because sup_neighbour itself may become negative if going over the 4-byte integer limit;
+!! DK DK but this should never happen in practice (by far)...
+ if( sup_neighbour > 2147483646 ) then
+ print *,'size exceeds integer 4-byte limit: ',sup_neighbour,nsize
+ print *,'bit size fortran: ',bit_size(sup_neighbour)
+ stop 'ERROR: sup_neighbour is too large'
endif
- sup_neighbour = NGNOD * nsize - (NGNOD + (NGNOD/2 - 1)*nfaces)
+ print *, ' nsize = ',nsize, 'sup_neighbour = ', sup_neighbour
- print*, ' nsize = ',nsize, 'sup_neighbour = ', sup_neighbour
-
end subroutine check_valence
!----------------------------------------------------------------------------------------------
@@ -639,7 +641,7 @@
implicit none
! local parameters
- integer, dimension(:),allocatable :: num_material
+ integer, dimension(:), allocatable :: num_material
integer :: ier
! starts from 0
@@ -658,9 +660,12 @@
call mesh2dual_ncommonnodes(nspec, nnodes, nsize, sup_neighbour, elmnts, xadj, adjncy, nnodes_elmnts, &
nodes_elmnts, max_neighbour, 1)
- print*, 'mesh2dual: '
+ print*, 'mesh2dual:'
print*, ' max_neighbour = ',max_neighbour
+!! DK DK Oct 2012: added this safety test
+ if(max_neighbour > sup_neighbour) stop 'found max_neighbour > sup_neighbour in domain decomposition'
+
nb_edges = xadj(nspec+1)
! allocates & initializes partioning of elements
@@ -683,12 +688,12 @@
! (which are counted then as elastic elements)
num_material(:) = mat(1,:)
- ! in case of acoustic/elastic/poro simulation, weights elements accordingly
+ ! in case of acoustic/elastic/poro simulation, assign different weights to elements accordingly
+!! DK DK Oct 2012: this should include CPML weights as well in the future
call acoustic_elastic_poro_load(elmnts_load,nspec,count_def_mat,count_undef_mat, &
num_material,mat_prop,undef_mat_prop)
-
! SCOTCH partitioning
! we use default strategy for partitioning, thus omit specifing explicit strategy .
@@ -776,7 +781,6 @@
! re-partitioning puts poroelastic-elastic coupled elements into same partition
! integer :: nfaces_coupled
! integer, dimension(:,:), pointer :: faces_coupled
-
! TODO: supposed to rebalance, but currently broken
call poro_elastic_repartitioning (nspec, nnodes, elmnts, &
count_def_mat, num_material , mat_prop, &
@@ -804,8 +808,6 @@
tab_size_interfaces, ninterfaces, &
nparts)
-
-
!or: uncomment if you want acoustic/elastic boundaries NOT to be separated into different MPI partitions
!call build_interfaces_no_ac_el_sep(nspec, sup_neighbour, part, elmnts, &
! xadj, adjncy, tab_interfaces, &
@@ -909,7 +911,7 @@
close(IIN_database)
- end do
+ enddo
print*, 'partitions: '
print*, ' num = ',nparts
print*
Modified: seismo/3D/SPECFEM3D/trunk/src/decompose_mesh_SCOTCH/part_decompose_mesh_SCOTCH.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/decompose_mesh_SCOTCH/part_decompose_mesh_SCOTCH.f90 2012-10-13 18:42:30 UTC (rev 20836)
+++ seismo/3D/SPECFEM3D/trunk/src/decompose_mesh_SCOTCH/part_decompose_mesh_SCOTCH.f90 2012-10-13 22:18:10 UTC (rev 20837)
@@ -30,10 +30,10 @@
include "../shared/constants.h"
-! Useful kind types
- integer ,parameter :: short = SELECTED_INT_KIND(4), long = SELECTED_INT_KIND(18)
+! useful kind types for short integer (4 bytes) and long integers (8 bytes)
+ integer ,parameter :: short = 4, long = 8
-! Number of faces per element.
+! number of faces per element.
integer, parameter :: nfaces = 6
! acoustic-elastic-poroelastic load balancing:
@@ -95,7 +95,7 @@
connectivity = 0
elem_base = nodes_elmnts(k+j*nsize)
elem_target = nodes_elmnts(l+j*nsize)
- do n = 1, NGNOD
+ do n = 1, NGNOD_EIGHT_CORNERS
num_node = elmnts(NGNOD*elem_base+n-1)
do m = 0, nnodes_elmnts(num_node)-1
if ( nodes_elmnts(m+num_node*nsize) == elem_target ) then
@@ -121,14 +121,14 @@
xadj(nodes_elmnts(k+j*nsize)) = xadj(nodes_elmnts(k+j*nsize)) + 1
if (xadj(nodes_elmnts(k+j*nsize))>sup_neighbour) &
- stop 'ERROR: too many neighbours per element, modify the mesh.'
+ stop 'ERROR: too many neighbours per element, error in the mesh or in the code.'
adjncy(nodes_elmnts(l+j*nsize)*sup_neighbour &
+ xadj(nodes_elmnts(l+j*nsize))) = nodes_elmnts(k+j*nsize)
xadj(nodes_elmnts(l+j*nsize)) = xadj(nodes_elmnts(l+j*nsize)) + 1
if (xadj(nodes_elmnts(l+j*nsize))>sup_neighbour) &
- stop 'ERROR: too many neighbours per element, modify the mesh.'
+ stop 'ERROR: too many neighbours per element, error in the mesh or in the code.'
end if
end if
end do
@@ -357,8 +357,8 @@
tab_interfaces(tab_size_interfaces(num_interface)*7+num_edge*7+0) = el
tab_interfaces(tab_size_interfaces(num_interface)*7+num_edge*7+1) = adjncy(el_adj)
ncommon_nodes = 0
- do num_node = 0, NGNOD-1
- do num_node_bis = 0, NGNOD-1
+ do num_node = 0, NGNOD_EIGHT_CORNERS-1
+ do num_node_bis = 0, NGNOD_EIGHT_CORNERS-1
if ( elmnts(el*NGNOD+num_node) == elmnts(adjncy(el_adj)*NGNOD+num_node_bis) ) then
tab_interfaces(tab_size_interfaces(num_interface)*7 &
+num_edge*7+3+ncommon_nodes) &
@@ -518,8 +518,8 @@
tab_interfaces(tab_size_interfaces(num_interface)*7+num_edge*7+0) = el
tab_interfaces(tab_size_interfaces(num_interface)*7+num_edge*7+1) = adjncy(el_adj)
ncommon_nodes = 0
- do num_node = 0, NGNOD-1
- do num_node_bis = 0, NGNOD-1
+ do num_node = 0, NGNOD_EIGHT_CORNERS-1
+ do num_node_bis = 0, NGNOD_EIGHT_CORNERS-1
if ( elmnts(el*NGNOD+num_node) == elmnts(adjncy(el_adj)*NGNOD+num_node_bis) ) then
tab_interfaces(tab_size_interfaces(num_interface)*7 &
+num_edge*7+3+ncommon_nodes) &
@@ -963,6 +963,8 @@
! format:
! # ispec_local # material_index_1 # material_index_2 # corner_id1 # corner_id2 # ... # corner_id8
+ ! or
+ ! # ispec_local # material_index_1 # material_index_2 # corner_id1 # corner_id2 # ... # corner_id27
write(IIN_database) glob2loc_elmnts(i)+1, num_modele(1,i+1), &
num_modele(2,i+1),(loc_nodes(k)+1, k=0,NGNOD-1)
end if
@@ -1451,8 +1453,7 @@
! partition index on each element
integer, dimension(0:nspec-1) :: part
- ! mesh element indexing
- ! ( elmnts(NGNOD,nspec) )
+ ! mesh element indexing ( elmnts(NGNOD,nspec) )
integer, dimension(0:NGNOD*nspec-1) :: elmnts
! moho surface
@@ -1509,14 +1510,14 @@
! loops over all element corners
counter = 0
- do i=0,NGNOD-1
+ do i=0,NGNOD_EIGHT_CORNERS-1
! note: assumes that node indices in elmnts array are in the range from 0 to nodes-1
inode = elmnts(el*NGNOD+i)
if( node_is_moho(inode) ) counter = counter + 1
enddo
! sets flag if it has a surface
- if( counter == 4 ) is_moho(el) = .true.
+ if( counter == NGNOD2D_FOUR_CORNERS ) is_moho(el) = .true.
enddo
! statistics output
Modified: seismo/3D/SPECFEM3D/trunk/src/generate_databases/create_regions_mesh.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/generate_databases/create_regions_mesh.f90 2012-10-13 18:42:30 UTC (rev 20836)
+++ seismo/3D/SPECFEM3D/trunk/src/generate_databases/create_regions_mesh.f90 2012-10-13 22:18:10 UTC (rev 20837)
@@ -272,11 +272,11 @@
! local parameters
integer :: ier
- allocate( xelm(NGNOD),yelm(NGNOD),zelm(NGNOD),stat=ier)
+ allocate(xelm(NGNOD),yelm(NGNOD),zelm(NGNOD),stat=ier)
if( ier /= 0 ) stop 'error allocating array xelm etc.'
- allocate( qmu_attenuation_store(NGLLX,NGLLY,NGLLZ,nspec),stat=ier)
- if(ier /= 0) call exit_MPI(myrank,'not enough memory to allocate arrays')
+ allocate(qmu_attenuation_store(NGLLX,NGLLY,NGLLZ,nspec),stat=ier)
+ if( ier /= 0 ) call exit_MPI(myrank,'not enough memory to allocate arrays')
! create the name for the database of the current slide and region
call create_name_database(prname,myrank,LOCAL_PATH)
@@ -291,63 +291,63 @@
! 3D shape functions and their derivatives
allocate(shape3D(NGNOD,NGLLX,NGLLY,NGLLZ), &
- dershape3D(NDIM,NGNOD,NGLLX,NGLLY,NGLLZ),stat=ier)
+ dershape3D(NDIM,NGNOD,NGLLX,NGLLY,NGLLZ),stat=ier)
if( ier /= 0 ) stop 'error allocating array shape3D etc.'
! 2D shape functions and their derivatives
allocate(shape2D_x(NGNOD2D,NGLLY,NGLLZ), &
- shape2D_y(NGNOD2D,NGLLX,NGLLZ), &
- shape2D_bottom(NGNOD2D,NGLLX,NGLLY), &
- shape2D_top(NGNOD2D,NGLLX,NGLLY),stat=ier)
+ shape2D_y(NGNOD2D,NGLLX,NGLLZ), &
+ shape2D_bottom(NGNOD2D,NGLLX,NGLLY), &
+ shape2D_top(NGNOD2D,NGLLX,NGLLY),stat=ier)
if( ier /= 0 ) stop 'error allocating array shape2D_x etc.'
allocate(dershape2D_x(NDIM2D,NGNOD2D,NGLLY,NGLLZ), &
- dershape2D_y(NDIM2D,NGNOD2D,NGLLX,NGLLZ), &
- dershape2D_bottom(NDIM2D,NGNOD2D,NGLLX,NGLLY), &
- dershape2D_top(NDIM2D,NGNOD2D,NGLLX,NGLLY),stat=ier)
+ dershape2D_y(NDIM2D,NGNOD2D,NGLLX,NGLLZ), &
+ dershape2D_bottom(NDIM2D,NGNOD2D,NGLLX,NGLLY), &
+ dershape2D_top(NDIM2D,NGNOD2D,NGLLX,NGLLY),stat=ier)
if( ier /= 0 ) stop 'error allocating array dershape2D_x etc.'
allocate(wgllwgll_xy(NGLLX,NGLLY), &
- wgllwgll_xz(NGLLX,NGLLZ), &
- wgllwgll_yz(NGLLY,NGLLZ),stat=ier)
+ wgllwgll_xz(NGLLX,NGLLZ), &
+ wgllwgll_yz(NGLLY,NGLLZ),stat=ier)
if(ier /= 0) call exit_MPI(myrank,'not enough memory to allocate arrays')
! Stacey
allocate(rho_vp(NGLLX,NGLLY,NGLLZ,nspec), &
- rho_vs(NGLLX,NGLLY,NGLLZ,nspec),stat=ier)
+ rho_vs(NGLLX,NGLLY,NGLLZ,nspec),stat=ier)
if(ier /= 0) call exit_MPI(myrank,'not enough memory to allocate arrays')
! array with model density
allocate(rhostore(NGLLX,NGLLY,NGLLZ,nspec), &
- kappastore(NGLLX,NGLLY,NGLLZ,nspec), &
- mustore(NGLLX,NGLLY,NGLLZ,nspec),stat=ier)
+ kappastore(NGLLX,NGLLY,NGLLZ,nspec), &
+ mustore(NGLLX,NGLLY,NGLLZ,nspec),stat=ier)
!vpstore(NGLLX,NGLLY,NGLLZ,nspec), &
!vsstore(NGLLX,NGLLY,NGLLZ,nspec),
if(ier /= 0) call exit_MPI(myrank,'not enough memory to allocate arrays')
! array with poroelastic model
allocate(rhoarraystore(2,NGLLX,NGLLY,NGLLZ,nspec), &
- kappaarraystore(3,NGLLX,NGLLY,NGLLZ,nspec), &
- etastore(NGLLX,NGLLY,NGLLZ,nspec), &
- tortstore(NGLLX,NGLLY,NGLLZ,nspec), &
- phistore(NGLLX,NGLLY,NGLLZ,nspec), &
- rho_vpI(NGLLX,NGLLY,NGLLZ,nspec), &
- rho_vpII(NGLLX,NGLLY,NGLLZ,nspec), &
- rho_vsI(NGLLX,NGLLY,NGLLZ,nspec), &
- permstore(6,NGLLX,NGLLY,NGLLZ,nspec), stat=ier)
+ kappaarraystore(3,NGLLX,NGLLY,NGLLZ,nspec), &
+ etastore(NGLLX,NGLLY,NGLLZ,nspec), &
+ tortstore(NGLLX,NGLLY,NGLLZ,nspec), &
+ phistore(NGLLX,NGLLY,NGLLZ,nspec), &
+ rho_vpI(NGLLX,NGLLY,NGLLZ,nspec), &
+ rho_vpII(NGLLX,NGLLY,NGLLZ,nspec), &
+ rho_vsI(NGLLX,NGLLY,NGLLZ,nspec), &
+ permstore(6,NGLLX,NGLLY,NGLLZ,nspec), stat=ier)
if(ier /= 0) call exit_MPI(myrank,'not enough memory to allocate arrays')
! arrays with mesh parameters
allocate(xixstore(NGLLX,NGLLY,NGLLZ,nspec), &
- xiystore(NGLLX,NGLLY,NGLLZ,nspec), &
- xizstore(NGLLX,NGLLY,NGLLZ,nspec), &
- etaxstore(NGLLX,NGLLY,NGLLZ,nspec), &
- etaystore(NGLLX,NGLLY,NGLLZ,nspec), &
- etazstore(NGLLX,NGLLY,NGLLZ,nspec), &
- gammaxstore(NGLLX,NGLLY,NGLLZ,nspec), &
- gammaystore(NGLLX,NGLLY,NGLLZ,nspec), &
- gammazstore(NGLLX,NGLLY,NGLLZ,nspec), &
- jacobianstore(NGLLX,NGLLY,NGLLZ,nspec),stat=ier)
+ xiystore(NGLLX,NGLLY,NGLLZ,nspec), &
+ xizstore(NGLLX,NGLLY,NGLLZ,nspec), &
+ etaxstore(NGLLX,NGLLY,NGLLZ,nspec), &
+ etaystore(NGLLX,NGLLY,NGLLZ,nspec), &
+ etazstore(NGLLX,NGLLY,NGLLZ,nspec), &
+ gammaxstore(NGLLX,NGLLY,NGLLZ,nspec), &
+ gammaystore(NGLLX,NGLLY,NGLLZ,nspec), &
+ gammazstore(NGLLX,NGLLY,NGLLZ,nspec), &
+ jacobianstore(NGLLX,NGLLY,NGLLZ,nspec),stat=ier)
if(ier /= 0) call exit_MPI(myrank,'not enough memory to allocate arrays')
! absorbing boundary
@@ -358,9 +358,9 @@
! allocates arrays to store info for each face (assumes NGLLX=NGLLY=NGLLZ)
allocate( abs_boundary_ispec(num_abs_boundary_faces), &
- abs_boundary_ijk(3,NGLLSQUARE,num_abs_boundary_faces), &
- abs_boundary_jacobian2Dw(NGLLSQUARE,num_abs_boundary_faces), &
- abs_boundary_normal(NDIM,NGLLSQUARE,num_abs_boundary_faces),stat=ier)
+ abs_boundary_ijk(3,NGLLSQUARE,num_abs_boundary_faces), &
+ abs_boundary_jacobian2Dw(NGLLSQUARE,num_abs_boundary_faces), &
+ abs_boundary_normal(NDIM,NGLLSQUARE,num_abs_boundary_faces),stat=ier)
if(ier /= 0) call exit_MPI(myrank,'not enough memory to allocate arrays')
! free surface faces
@@ -368,9 +368,9 @@
! allocates arrays to store info for each face (assumes NGLLX=NGLLY=NGLLZ)
allocate( free_surface_ispec(num_free_surface_faces), &
- free_surface_ijk(3,NGLLSQUARE,num_free_surface_faces), &
- free_surface_jacobian2Dw(NGLLSQUARE,num_free_surface_faces), &
- free_surface_normal(NDIM,NGLLSQUARE,num_free_surface_faces),stat=ier)
+ free_surface_ijk(3,NGLLSQUARE,num_free_surface_faces), &
+ free_surface_jacobian2Dw(NGLLSQUARE,num_free_surface_faces), &
+ free_surface_normal(NDIM,NGLLSQUARE,num_free_surface_faces),stat=ier)
if(ier /= 0) call exit_MPI(myrank,'not enough memory to allocate arrays')
! array with anisotropy
@@ -380,32 +380,32 @@
NSPEC_ANISO = 1
endif
allocate(c11store(NGLLX,NGLLY,NGLLZ,NSPEC_ANISO), &
- c12store(NGLLX,NGLLY,NGLLZ,NSPEC_ANISO), &
- c13store(NGLLX,NGLLY,NGLLZ,NSPEC_ANISO), &
- c14store(NGLLX,NGLLY,NGLLZ,NSPEC_ANISO), &
- c15store(NGLLX,NGLLY,NGLLZ,NSPEC_ANISO), &
- c16store(NGLLX,NGLLY,NGLLZ,NSPEC_ANISO), &
- c22store(NGLLX,NGLLY,NGLLZ,NSPEC_ANISO), &
- c23store(NGLLX,NGLLY,NGLLZ,NSPEC_ANISO), &
- c24store(NGLLX,NGLLY,NGLLZ,NSPEC_ANISO), &
- c25store(NGLLX,NGLLY,NGLLZ,NSPEC_ANISO), &
- c26store(NGLLX,NGLLY,NGLLZ,NSPEC_ANISO), &
- c33store(NGLLX,NGLLY,NGLLZ,NSPEC_ANISO), &
- c34store(NGLLX,NGLLY,NGLLZ,NSPEC_ANISO), &
- c35store(NGLLX,NGLLY,NGLLZ,NSPEC_ANISO), &
- c36store(NGLLX,NGLLY,NGLLZ,NSPEC_ANISO), &
- c44store(NGLLX,NGLLY,NGLLZ,NSPEC_ANISO), &
- c45store(NGLLX,NGLLY,NGLLZ,NSPEC_ANISO), &
- c46store(NGLLX,NGLLY,NGLLZ,NSPEC_ANISO), &
- c55store(NGLLX,NGLLY,NGLLZ,NSPEC_ANISO), &
- c56store(NGLLX,NGLLY,NGLLZ,NSPEC_ANISO), &
- c66store(NGLLX,NGLLY,NGLLZ,NSPEC_ANISO),stat=ier)
+ c12store(NGLLX,NGLLY,NGLLZ,NSPEC_ANISO), &
+ c13store(NGLLX,NGLLY,NGLLZ,NSPEC_ANISO), &
+ c14store(NGLLX,NGLLY,NGLLZ,NSPEC_ANISO), &
+ c15store(NGLLX,NGLLY,NGLLZ,NSPEC_ANISO), &
+ c16store(NGLLX,NGLLY,NGLLZ,NSPEC_ANISO), &
+ c22store(NGLLX,NGLLY,NGLLZ,NSPEC_ANISO), &
+ c23store(NGLLX,NGLLY,NGLLZ,NSPEC_ANISO), &
+ c24store(NGLLX,NGLLY,NGLLZ,NSPEC_ANISO), &
+ c25store(NGLLX,NGLLY,NGLLZ,NSPEC_ANISO), &
+ c26store(NGLLX,NGLLY,NGLLZ,NSPEC_ANISO), &
+ c33store(NGLLX,NGLLY,NGLLZ,NSPEC_ANISO), &
+ c34store(NGLLX,NGLLY,NGLLZ,NSPEC_ANISO), &
+ c35store(NGLLX,NGLLY,NGLLZ,NSPEC_ANISO), &
+ c36store(NGLLX,NGLLY,NGLLZ,NSPEC_ANISO), &
+ c44store(NGLLX,NGLLY,NGLLZ,NSPEC_ANISO), &
+ c45store(NGLLX,NGLLY,NGLLZ,NSPEC_ANISO), &
+ c46store(NGLLX,NGLLY,NGLLZ,NSPEC_ANISO), &
+ c55store(NGLLX,NGLLY,NGLLZ,NSPEC_ANISO), &
+ c56store(NGLLX,NGLLY,NGLLZ,NSPEC_ANISO), &
+ c66store(NGLLX,NGLLY,NGLLZ,NSPEC_ANISO),stat=ier)
if(ier /= 0) call exit_MPI(myrank,'not enough memory to allocate arrays')
! material flags
allocate( ispec_is_acoustic(nspec), &
- ispec_is_elastic(nspec), &
- ispec_is_poroelastic(nspec), stat=ier)
+ ispec_is_elastic(nspec), &
+ ispec_is_poroelastic(nspec), stat=ier)
if(ier /= 0) call exit_MPI(myrank,'not enough memory to allocate arrays')
! initializes Moho surface
@@ -626,7 +626,7 @@
integer, dimension(NGLLX,NGLLY,NGLLZ,nspec) :: ibool
! local parameters
- real(kind=CUSTOM_REAL),dimension(NGNOD2D) :: xcoord,ycoord,zcoord
+ real(kind=CUSTOM_REAL),dimension(NGNOD2D_FOUR_CORNERS) :: xcoord,ycoord,zcoord
real(kind=CUSTOM_REAL) :: jacobian2Dw_face(NGLLX,NGLLY)
real(kind=CUSTOM_REAL) :: normal_face(NDIM,NGLLX,NGLLY)
real(kind=CUSTOM_REAL),dimension(NDIM):: normal
@@ -679,7 +679,7 @@
! determines element face by given CUBIT corners
! (note: uses point locations rather than point indices to find the element face,
! because the indices refer no more to the newly indexed ibool array )
- do icorner=1,NGNOD2D
+ do icorner=1,NGNOD2D_FOUR_CORNERS
xcoord(icorner) = nodes_coords_ext_mesh(1,nodes_ibelm_moho(icorner,ispec2D))
ycoord(icorner) = nodes_coords_ext_mesh(2,nodes_ibelm_moho(icorner,ispec2D))
zcoord(icorner) = nodes_coords_ext_mesh(3,nodes_ibelm_moho(icorner,ispec2D))
@@ -699,7 +699,7 @@
xstore_dummy,ystore_dummy,zstore_dummy,ibool,nglob_dummy,&
dershape2D_x,dershape2D_y,dershape2D_bottom,dershape2D_top, &
wgllwgll_xy,wgllwgll_xz,wgllwgll_yz,&
- ispec,iface,jacobian2Dw_face,normal_face,NGLLX,NGLLZ)
+ ispec,iface,jacobian2Dw_face,normal_face,NGLLX,NGLLZ,NGNOD2D)
! normal convention: points away from element
! switch normal direction if necessary
@@ -755,7 +755,7 @@
do iface = 1,6
! checks if corners of face on surface
counter = 0
- do icorner = 1,NGNOD2D
+ do icorner = 1,NGNOD2D_FOUR_CORNERS
i = iface_all_corner_ijk(1,icorner,iface)
j = iface_all_corner_ijk(2,icorner,iface)
k = iface_all_corner_ijk(3,icorner,iface)
@@ -773,7 +773,7 @@
enddo
! stores moho informations
- if( counter == NGNOD2D ) then
+ if( counter == NGNOD2D_FOUR_CORNERS ) then
! gets face GLL points i,j,k indices from element face
call get_element_face_gll_indices(iface,ijk_face,NGLLX,NGLLY)
@@ -784,7 +784,7 @@
xstore_dummy,ystore_dummy,zstore_dummy,ibool,nglob_dummy,&
dershape2D_x,dershape2D_y,dershape2D_bottom,dershape2D_top, &
wgllwgll_xy,wgllwgll_xz,wgllwgll_yz,&
- ispec,iface,jacobian2Dw_face,normal_face,NGLLX,NGLLZ)
+ ispec,iface,jacobian2Dw_face,normal_face,NGLLX,NGLLZ,NGNOD2D)
! normal convention: points away from element
! switch normal direction if necessary
Modified: seismo/3D/SPECFEM3D/trunk/src/generate_databases/generate_databases.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/generate_databases/generate_databases.f90 2012-10-13 18:42:30 UTC (rev 20836)
+++ seismo/3D/SPECFEM3D/trunk/src/generate_databases/generate_databases.f90 2012-10-13 22:18:10 UTC (rev 20837)
@@ -333,9 +333,12 @@
if(CUSTOM_REAL /= SIZE_REAL .and. CUSTOM_REAL /= SIZE_DOUBLE) &
call exit_MPI(myrank,'wrong size of CUSTOM_REAL for reals')
- 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')
+ if(NGNOD /= 8 .and. NGNOD /= 27) call exit_MPI(myrank,'number of element control nodes must be 8 or 27')
+ if(NGNOD2D /= 4 .and. NGNOD2D /= 9) call exit_MPI(myrank,'number of face control nodes must be 4 or 9')
+ if(NGNOD == 8 .and. NGNOD2D /= 4) call exit_MPI(myrank,'elements with 8 control nodes must have NGNOD2D == 4')
+ if(NGNOD == 27 .and. NGNOD2D /= 9) call exit_MPI(myrank,'elements with 27 control nodes must have NGNOD2D == 9')
+
! 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/generate_databases/get_absorbing_boundary.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/generate_databases/get_absorbing_boundary.f90 2012-10-13 18:42:30 UTC (rev 20836)
+++ seismo/3D/SPECFEM3D/trunk/src/generate_databases/get_absorbing_boundary.f90 2012-10-13 22:18:10 UTC (rev 20837)
@@ -111,7 +111,7 @@
xstore_dummy,ystore_dummy,zstore_dummy,ibool,nglob_dummy,&
dershape2D_x,dershape2D_y,dershape2D_bottom,dershape2D_top, &
wgllwgll_xy,wgllwgll_xz,wgllwgll_yz,&
- ispec,iface,jacobian2Dw_face,normal_face,NGLLX,NGLLZ)
+ ispec,iface,jacobian2Dw_face,normal_face,NGLLX,NGLLZ,NGNOD2D)
! normal convention: points away from element
! switch normal direction if necessary
@@ -173,7 +173,7 @@
xstore_dummy,ystore_dummy,zstore_dummy,ibool,nglob_dummy,&
dershape2D_x,dershape2D_y,dershape2D_bottom,dershape2D_top, &
wgllwgll_xy,wgllwgll_xz,wgllwgll_yz,&
- ispec,iface,jacobian2Dw_face,normal_face,NGLLX,NGLLZ)
+ ispec,iface,jacobian2Dw_face,normal_face,NGLLX,NGLLZ,NGNOD2D)
! normal convention: points away from element
! switch normal direction if necessary
@@ -235,7 +235,7 @@
xstore_dummy,ystore_dummy,zstore_dummy,ibool,nglob_dummy,&
dershape2D_x,dershape2D_y,dershape2D_bottom,dershape2D_top, &
wgllwgll_xy,wgllwgll_xz,wgllwgll_yz,&
- ispec,iface,jacobian2Dw_face,normal_face,NGLLY,NGLLZ)
+ ispec,iface,jacobian2Dw_face,normal_face,NGLLY,NGLLZ,NGNOD2D)
! normal convention: points away from element
! switch normal direction if necessary
@@ -297,7 +297,7 @@
xstore_dummy,ystore_dummy,zstore_dummy,ibool,nglob_dummy, &
dershape2D_x,dershape2D_y,dershape2D_bottom,dershape2D_top, &
wgllwgll_xy,wgllwgll_xz,wgllwgll_yz, &
- ispec,iface,jacobian2Dw_face,normal_face,NGLLY,NGLLZ)
+ ispec,iface,jacobian2Dw_face,normal_face,NGLLY,NGLLZ,NGNOD2D)
! normal convention: points away from element
! switch normal direction if necessary
@@ -359,7 +359,7 @@
xstore_dummy,ystore_dummy,zstore_dummy,ibool,nglob_dummy, &
dershape2D_x,dershape2D_y,dershape2D_bottom,dershape2D_top, &
wgllwgll_xy,wgllwgll_xz,wgllwgll_yz, &
- ispec,iface,jacobian2Dw_face,normal_face,NGLLX,NGLLY)
+ ispec,iface,jacobian2Dw_face,normal_face,NGLLX,NGLLY,NGNOD2D)
! normal convention: points away from element
! switch normal direction if necessary
@@ -423,7 +423,7 @@
xstore_dummy,ystore_dummy,zstore_dummy,ibool,nglob_dummy, &
dershape2D_x,dershape2D_y,dershape2D_bottom,dershape2D_top, &
wgllwgll_xy,wgllwgll_xz,wgllwgll_yz, &
- ispec,iface,jacobian2Dw_face,normal_face,NGLLX,NGLLY)
+ ispec,iface,jacobian2Dw_face,normal_face,NGLLX,NGLLY,NGNOD2D)
! normal convention: points away from element
! switch normal direction if necessary
Modified: seismo/3D/SPECFEM3D/trunk/src/generate_databases/get_coupling_surfaces.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/generate_databases/get_coupling_surfaces.f90 2012-10-13 18:42:30 UTC (rev 20836)
+++ seismo/3D/SPECFEM3D/trunk/src/generate_databases/get_coupling_surfaces.f90 2012-10-13 22:18:10 UTC (rev 20837)
@@ -289,7 +289,7 @@
xstore_dummy,ystore_dummy,zstore_dummy,ibool,nglob_dummy, &
dershape2D_x,dershape2D_y,dershape2D_bottom,dershape2D_top, &
wgllwgll_xy,wgllwgll_xz,wgllwgll_yz, &
- ispec,iface_ref,jacobian2Dw_face,normal_face,NGLLX,NGLLY)
+ ispec,iface_ref,jacobian2Dw_face,normal_face,NGLLX,NGLLY,NGNOD2D)
! normal convention: points away from acoustic, reference element
! switch normal direction if necessary
@@ -471,7 +471,7 @@
xstore_dummy,ystore_dummy,zstore_dummy,ibool,nglob_dummy, &
dershape2D_x,dershape2D_y,dershape2D_bottom,dershape2D_top, &
wgllwgll_xy,wgllwgll_xz,wgllwgll_yz, &
- ispec,iface_ref,jacobian2Dw_face,normal_face,NGLLX,NGLLY)
+ ispec,iface_ref,jacobian2Dw_face,normal_face,NGLLX,NGLLY,NGNOD2D)
! normal convention: points away from acoustic, reference element
! switch normal direction if necessary
@@ -668,7 +668,7 @@
xstore_dummy,ystore_dummy,zstore_dummy,ibool,nglob_dummy, &
dershape2D_x,dershape2D_y,dershape2D_bottom,dershape2D_top, &
wgllwgll_xy,wgllwgll_xz,wgllwgll_yz, &
- ispec,iface_ref,jacobian2Dw_face,normal_face,NGLLX,NGLLY)
+ ispec,iface_ref,jacobian2Dw_face,normal_face,NGLLX,NGLLY,NGNOD2D)
! normal convention: points away from poroelastic, reference element
do j=1,NGLLY
Modified: seismo/3D/SPECFEM3D/trunk/src/generate_databases/read_partition_files.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/generate_databases/read_partition_files.f90 2012-10-13 18:42:30 UTC (rev 20836)
+++ seismo/3D/SPECFEM3D/trunk/src/generate_databases/read_partition_files.f90 2012-10-13 22:18:10 UTC (rev 20837)
@@ -115,9 +115,7 @@
! format:
! # ispec_local # material_index_1 # material_index_2 # corner_id1 # corner_id2 # ... # corner_id8
- read(IIN) dummy_elmnt, mat_ext_mesh(1,ispec),mat_ext_mesh(2,ispec), &
- elmnts_ext_mesh(1,ispec), elmnts_ext_mesh(2,ispec), elmnts_ext_mesh(3,ispec), elmnts_ext_mesh(4,ispec), &
- elmnts_ext_mesh(5,ispec), elmnts_ext_mesh(6,ispec), elmnts_ext_mesh(7,ispec), elmnts_ext_mesh(8,ispec)
+ read(IIN) dummy_elmnt, mat_ext_mesh(1,ispec),mat_ext_mesh(2,ispec),(elmnts_ext_mesh(j,ispec),j=1,NGNOD)
! check debug
if( dummy_elmnt /= ispec) stop "error ispec order in materials file"
@@ -153,40 +151,40 @@
NSPEC2D_BOTTOM = nspec2D_bottom_ext
NSPEC2D_TOP = nspec2D_top_ext
- allocate(ibelm_xmin(nspec2D_xmin),nodes_ibelm_xmin(4,nspec2D_xmin),stat=ier)
+ allocate(ibelm_xmin(nspec2D_xmin),nodes_ibelm_xmin(NGNOD2D,nspec2D_xmin),stat=ier)
if( ier /= 0 ) stop 'error allocating array ibelm_xmin etc.'
do ispec2D = 1,nspec2D_xmin
- read(IIN) ibelm_xmin(ispec2D),(nodes_ibelm_xmin(j,ispec2D),j=1,4)
+ read(IIN) ibelm_xmin(ispec2D),(nodes_ibelm_xmin(j,ispec2D),j=1,NGNOD2D)
end do
- allocate(ibelm_xmax(nspec2D_xmax),nodes_ibelm_xmax(4,nspec2D_xmax),stat=ier)
+ allocate(ibelm_xmax(nspec2D_xmax),nodes_ibelm_xmax(NGNOD2D,nspec2D_xmax),stat=ier)
if( ier /= 0 ) stop 'error allocating array ibelm_xmax etc.'
do ispec2D = 1,nspec2D_xmax
- read(IIN) ibelm_xmax(ispec2D),(nodes_ibelm_xmax(j,ispec2D),j=1,4)
+ read(IIN) ibelm_xmax(ispec2D),(nodes_ibelm_xmax(j,ispec2D),j=1,NGNOD2D)
end do
- allocate(ibelm_ymin(nspec2D_ymin),nodes_ibelm_ymin(4,nspec2D_ymin),stat=ier)
+ allocate(ibelm_ymin(nspec2D_ymin),nodes_ibelm_ymin(NGNOD2D,nspec2D_ymin),stat=ier)
if( ier /= 0 ) stop 'error allocating array ibelm_ymin'
do ispec2D = 1,nspec2D_ymin
- read(IIN) ibelm_ymin(ispec2D),(nodes_ibelm_ymin(j,ispec2D),j=1,4)
+ read(IIN) ibelm_ymin(ispec2D),(nodes_ibelm_ymin(j,ispec2D),j=1,NGNOD2D)
end do
- allocate(ibelm_ymax(nspec2D_ymax),nodes_ibelm_ymax(4,nspec2D_ymax),stat=ier)
+ allocate(ibelm_ymax(nspec2D_ymax),nodes_ibelm_ymax(NGNOD2D,nspec2D_ymax),stat=ier)
if( ier /= 0 ) stop 'error allocating array ibelm_ymax etc.'
do ispec2D = 1,nspec2D_ymax
- read(IIN) ibelm_ymax(ispec2D),(nodes_ibelm_ymax(j,ispec2D),j=1,4)
+ read(IIN) ibelm_ymax(ispec2D),(nodes_ibelm_ymax(j,ispec2D),j=1,NGNOD2D)
end do
- allocate(ibelm_bottom(nspec2D_bottom_ext),nodes_ibelm_bottom(4,nspec2D_bottom_ext),stat=ier)
+ allocate(ibelm_bottom(nspec2D_bottom_ext),nodes_ibelm_bottom(NGNOD2D,nspec2D_bottom_ext),stat=ier)
if( ier /= 0 ) stop 'error allocating array ibelm_bottom etc.'
do ispec2D = 1,nspec2D_bottom_ext
- read(IIN) ibelm_bottom(ispec2D),(nodes_ibelm_bottom(j,ispec2D),j=1,4)
+ read(IIN) ibelm_bottom(ispec2D),(nodes_ibelm_bottom(j,ispec2D),j=1,NGNOD2D)
end do
- allocate(ibelm_top(nspec2D_top_ext),nodes_ibelm_top(4,nspec2D_top_ext),stat=ier)
+ allocate(ibelm_top(nspec2D_top_ext),nodes_ibelm_top(NGNOD2D,nspec2D_top_ext),stat=ier)
if( ier /= 0 ) stop 'error allocating array ibelm_top etc.'
do ispec2D = 1,nspec2D_top_ext
- read(IIN) ibelm_top(ispec2D),(nodes_ibelm_top(j,ispec2D),j=1,4)
+ read(IIN) ibelm_top(ispec2D),(nodes_ibelm_top(j,ispec2D),j=1,NGNOD2D)
end do
call sum_all_i(nspec2D_xmin,num_xmin)
@@ -269,11 +267,11 @@
if( num_moho == 0 ) call exit_mpi(myrank,'error no moho mesh in database')
! reads in element informations
- allocate(ibelm_moho(nspec2D_moho_ext),nodes_ibelm_moho(4,nspec2D_moho_ext),stat=ier)
+ allocate(ibelm_moho(nspec2D_moho_ext),nodes_ibelm_moho(NGNOD2D,nspec2D_moho_ext),stat=ier)
if( ier /= 0 ) stop 'error allocating array ibelm_moho etc.'
do ispec2D = 1,nspec2D_moho_ext
! format: #element_id #node_id1 #node_id2 #node_id3 #node_id4
- read(IIN) ibelm_moho(ispec2D),(nodes_ibelm_moho(j,ispec2D),j=1,4)
+ read(IIN) ibelm_moho(ispec2D),(nodes_ibelm_moho(j,ispec2D),j=1,NGNOD2D)
end do
! user output
@@ -284,7 +282,7 @@
else
! allocate dummy array
nspec2D_moho_ext = 0
- allocate(ibelm_moho(nspec2D_moho_ext),nodes_ibelm_moho(4,nspec2D_moho_ext),stat=ier)
+ allocate(ibelm_moho(nspec2D_moho_ext),nodes_ibelm_moho(NGNOD2D,nspec2D_moho_ext),stat=ier)
if( ier /= 0 ) stop 'error allocating dumy array ibelm_moho etc.'
endif
Modified: seismo/3D/SPECFEM3D/trunk/src/shared/get_element_face.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/shared/get_element_face.f90 2012-10-13 18:42:30 UTC (rev 20836)
+++ seismo/3D/SPECFEM3D/trunk/src/shared/get_element_face.f90 2012-10-13 22:18:10 UTC (rev 20837)
@@ -288,6 +288,14 @@
integer :: ispec,iface,nspec,nglob
! face corner locations
+!! DK DK Oct 2012: in principle we should use NGNOD2D instead of NGNOD2D_FOUR_CORNERS when
+!! DK DK Oct 2012: computing the normal in the case of HEX27 elements, to be more precise;
+!! DK DK Oct 2012: but the code below would be a bit difficult to modify therefore we keep
+!! DK DK Oct 2012: using NGNOD2D_FOUR_CORNERS only for now.
+!! DK DK Oct 2012: If the face is flat (for instance along an absorbing edge) then the above
+!! DK DK Oct 2012: code is still exact even for HEX27, but if there is bathymetry for instance
+!! DK DK Oct 2012: along a curved fluid-solid interface then the code below is not as precise
+!! DK DK Oct 2012: as it should for HEX27.
real(kind=CUSTOM_REAL),dimension(NGNOD2D_FOUR_CORNERS) :: xcoord,ycoord,zcoord
! index array
@@ -317,7 +325,7 @@
face_n(:) = face_n(:)/tmp
! checks that this normal direction is outwards of element:
- ! takes additional corner out of face plane and determines scalarproduct to normal
+ ! takes additional corner out of face plane and determines scalar product (dot product) to normal
iglob = 0
select case( iface )
case(1) ! opposite to xmin face
@@ -340,10 +348,10 @@
v_tmp(2) = ystore_dummy(iglob) - ycoord(1)
v_tmp(3) = zstore_dummy(iglob) - zcoord(1)
- ! scalarproduct
+ ! scalar product (dot product)
tmp = v_tmp(1)*face_n(1) + v_tmp(2)*face_n(2) + v_tmp(3)*face_n(3)
- ! makes sure normal points outwards, that is points away from this additional corner and scalarproduct is negative
+ ! makes sure normal points outwards, that is points away from this additional corner and scalar product (dot product) is negative
if( tmp > 0.0_CUSTOM_REAL ) then
face_n(:) = - face_n(:)
endif
@@ -351,7 +359,6 @@
! in case given normal has zero length, sets it to computed face normal
! note: to avoid floating-point exception we use dble()
! values of normal(:) could be very small, almost zero, and lead to underflow
- !tmp = normal(1)*normal(1) + normal(2)*normal(2) + normal(3)*normal(3)
tmp = sngl(dble(normal(1))**2 + dble(normal(2))**2 + dble(normal(3))**2)
if( tmp < TINYVAL ) then
normal(:) = face_n(:)
@@ -359,7 +366,6 @@
endif
! otherwise determines orientation of normal and flips direction such that normal points outwards
-
tmp = face_n(1)*normal(1) + face_n(2)*normal(2) + face_n(3)*normal(3)
if( tmp < 0.0_CUSTOM_REAL ) then
!swap
@@ -388,6 +394,14 @@
integer :: ispec,iface,nspec,nglob
! face corner locations
+!! DK DK Oct 2012: in principle we should use NGNOD2D instead of NGNOD2D_FOUR_CORNERS when
+!! DK DK Oct 2012: computing the normal in the case of HEX27 elements, to be more precise;
+!! DK DK Oct 2012: but the code below would be a bit difficult to modify therefore we keep
+!! DK DK Oct 2012: using NGNOD2D_FOUR_CORNERS only for now.
+!! DK DK Oct 2012: If the face is flat (for instance along an absorbing edge) then the above
+!! DK DK Oct 2012: code is still exact even for HEX27, but if there is bathymetry for instance
+!! DK DK Oct 2012: along a curved fluid-solid interface then the code below is not as precise
+!! DK DK Oct 2012: as it should for HEX27.
real(kind=CUSTOM_REAL),dimension(NGNOD2D_FOUR_CORNERS) :: xcoord,ycoord,zcoord
! index array
@@ -420,7 +434,7 @@
face_n(:) = face_n(:)/tmp
! checks that this normal direction is outwards of element:
- ! takes additional corner out of face plane and determines scalarproduct to normal
+ ! takes additional corner out of face plane and determines scalar product (dot product) to normal
iglob = 0
select case( iface )
case(1) ! opposite to xmin face
@@ -443,10 +457,10 @@
v_tmp(2) = ystore_dummy(iglob) - ycoord(1)
v_tmp(3) = zstore_dummy(iglob) - zcoord(1)
- ! scalarproduct
+ ! scalar product (dot product)
tmp = v_tmp(1)*face_n(1) + v_tmp(2)*face_n(2) + v_tmp(3)*face_n(3)
- ! makes sure normal points outwards, that is points away from this additional corner and scalarproduct is negative
+ ! makes sure normal points outwards, that is points away from this additional corner and scalar product (dot product) is negative
if( tmp > 0.0 ) then
face_n(:) = - face_n(:)
endif
@@ -494,6 +508,14 @@
real(kind=CUSTOM_REAL) :: xstore_dummy(nglob),ystore_dummy(nglob),zstore_dummy(nglob)
! assumes NGNOD2D_FOUR_CORNERS == 4
+!! DK DK Oct 2012: in principle we should use NGNOD2D instead of NGNOD2D_FOUR_CORNERS when
+!! DK DK Oct 2012: computing the normal in the case of HEX27 elements, to be more precise;
+!! DK DK Oct 2012: but the code below would be a bit difficult to modify therefore we keep
+!! DK DK Oct 2012: using NGNOD2D_FOUR_CORNERS only for now.
+!! DK DK Oct 2012: If the face is flat (for instance along an absorbing edge) then the above
+!! DK DK Oct 2012: code is still exact even for HEX27, but if there is bathymetry for instance
+!! DK DK Oct 2012: along a curved fluid-solid interface then the code below is not as precise
+!! DK DK Oct 2012: as it should for HEX27.
integer,dimension(3,4,6) :: iface_all_corner_ijk
! local parameters
Modified: seismo/3D/SPECFEM3D/trunk/src/shared/get_force.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/shared/get_force.f90 2012-10-13 18:42:30 UTC (rev 20836)
+++ seismo/3D/SPECFEM3D/trunk/src/shared/get_force.f90 2012-10-13 22:18:10 UTC (rev 20837)
@@ -99,15 +99,15 @@
read(1,"(a)") string
read(string(21:len_trim(string)),*) factor_force_source(isource)
- ! read direction vector's East component
+ ! read direction vector's East component
read(1,"(a)") string
read(string(29:len_trim(string)),*) comp_dir_vect_source_E(isource)
- ! read direction vector's North component
+ ! read direction vector's North component
read(1,"(a)") string
read(string(29:len_trim(string)),*) comp_dir_vect_source_N(isource)
- ! read direction vector's vertical component
+ ! read direction vector's vertical component
read(1,"(a)") string
read(string(32:len_trim(string)),*) comp_dir_vect_source_Z_UP(isource)
@@ -125,7 +125,7 @@
endif
do isource=1,NSOURCES
-
+
! checks half-duration
! half-duration is the dominant frequency of the source
! point forces use a Ricker source time function
@@ -133,7 +133,7 @@
! (see constants.h: TINYVAL = 1.d-9 )
if( hdur(isource) < TINYVAL ) hdur(isource) = TINYVAL
- ! check (inclined) force source's direction vector
+ ! check (inclined) force source's direction vector
if( comp_dir_vect_source_E(isource) .eq. 0.d0 .and. comp_dir_vect_source_N(isource) .eq. 0.d0 &
.and. comp_dir_vect_source_Z_UP(isource) .eq. 0.d0 ) &
stop 'When using USE_FORCE_POINT_SOURCE make sure all forces have a non null direction vector'
Modified: seismo/3D/SPECFEM3D/trunk/src/shared/get_jacobian_boundaries.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/shared/get_jacobian_boundaries.f90 2012-10-13 18:42:30 UTC (rev 20836)
+++ seismo/3D/SPECFEM3D/trunk/src/shared/get_jacobian_boundaries.f90 2012-10-13 22:18:10 UTC (rev 20837)
@@ -28,7 +28,7 @@
xstore_dummy,ystore_dummy,zstore_dummy,ibool,nglob,&
dershape2D_x,dershape2D_y,dershape2D_bottom,dershape2D_top, &
wgllwgll_xy,wgllwgll_xz,wgllwgll_yz,&
- ispec,iface,jacobian2Dw_face,normal_face,NGLLA,NGLLB)
+ ispec,iface,jacobian2Dw_face,normal_face,NGLLA,NGLLB,NGNOD2D_value)
! returns jacobian2Dw_face and normal_face (pointing outwards of element)
@@ -47,10 +47,11 @@
real(kind=CUSTOM_REAL), dimension(NGLLA,NGLLB) :: jacobian2Dw_face
real(kind=CUSTOM_REAL), dimension(NDIM,NGLLA,NGLLB) :: normal_face
- 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)
+ integer :: NGNOD2D_value
+ double precision dershape2D_x(NDIM2D,NGNOD2D_value,NGLLY,NGLLZ)
+ double precision dershape2D_y(NDIM2D,NGNOD2D_value,NGLLX,NGLLZ)
+ double precision dershape2D_bottom(NDIM2D,NGNOD2D_value,NGLLX,NGLLY)
+ double precision dershape2D_top(NDIM2D,NGNOD2D_value,NGLLX,NGLLY)
double precision, dimension(NGLLX,NGLLY) :: wgllwgll_xy
double precision, dimension(NGLLX,NGLLZ) :: wgllwgll_xz
@@ -58,11 +59,10 @@
! local parameters
! face corners
- double precision xelm(NGNOD2D),yelm(NGNOD2D),zelm(NGNOD2D)
+ double precision xelm(NGNOD2D_value),yelm(NGNOD2D_value),zelm(NGNOD2D_value)
! 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')
+ if(NGNOD2D_value /= 4 .and. NGNOD2D_value /= 9) call exit_MPI(myrank,'surface elements should have 4 or 9 control nodes')
select case ( iface )
! on reference face: xmin
@@ -80,6 +80,24 @@
yelm(4)=ystore_dummy( ibool(1,1,NGLLZ,ispec) )
zelm(4)=zstore_dummy( ibool(1,1,NGLLZ,ispec) )
+ if(NGNOD2D_value == 9) then
+ xelm(5)=xstore_dummy( ibool(1,(NGLLY+1)/2,1,ispec) )
+ yelm(5)=ystore_dummy( ibool(1,(NGLLY+1)/2,1,ispec) )
+ zelm(5)=zstore_dummy( ibool(1,(NGLLY+1)/2,1,ispec) )
+ xelm(6)=xstore_dummy( ibool(1,NGLLY,(NGLLZ+1)/2,ispec) )
+ yelm(6)=ystore_dummy( ibool(1,NGLLY,(NGLLZ+1)/2,ispec) )
+ zelm(6)=zstore_dummy( ibool(1,NGLLY,(NGLLZ+1)/2,ispec) )
+ xelm(7)=xstore_dummy( ibool(1,(NGLLY+1)/2,NGLLZ,ispec) )
+ yelm(7)=ystore_dummy( ibool(1,(NGLLY+1)/2,NGLLZ,ispec) )
+ zelm(7)=zstore_dummy( ibool(1,(NGLLY+1)/2,NGLLZ,ispec) )
+ xelm(8)=xstore_dummy( ibool(1,1,(NGLLZ+1)/2,ispec) )
+ yelm(8)=ystore_dummy( ibool(1,1,(NGLLZ+1)/2,ispec) )
+ zelm(8)=zstore_dummy( ibool(1,1,(NGLLZ+1)/2,ispec) )
+ xelm(9)=xstore_dummy( ibool(1,(NGLLY+1)/2,(NGLLZ+1)/2,ispec) )
+ yelm(9)=ystore_dummy( ibool(1,(NGLLY+1)/2,(NGLLZ+1)/2,ispec) )
+ zelm(9)=zstore_dummy( ibool(1,(NGLLY+1)/2,(NGLLZ+1)/2,ispec) )
+ endif
+
call compute_jacobian_2D_face(myrank,xelm,yelm,zelm, &
dershape2D_x,wgllwgll_yz, &
jacobian2Dw_face,normal_face,NGLLY,NGLLZ)
@@ -99,6 +117,24 @@
yelm(4)=ystore_dummy( ibool(NGLLX,1,NGLLZ,ispec) )
zelm(4)=zstore_dummy( ibool(NGLLX,1,NGLLZ,ispec) )
+ if(NGNOD2D_value == 9) then
+ xelm(5)=xstore_dummy( ibool(NGLLX,(NGLLY+1)/2,1,ispec) )
+ yelm(5)=ystore_dummy( ibool(NGLLX,(NGLLY+1)/2,1,ispec) )
+ zelm(5)=zstore_dummy( ibool(NGLLX,(NGLLY+1)/2,1,ispec) )
+ xelm(6)=xstore_dummy( ibool(NGLLX,NGLLY,(NGLLZ+1)/2,ispec) )
+ yelm(6)=ystore_dummy( ibool(NGLLX,NGLLY,(NGLLZ+1)/2,ispec) )
+ zelm(6)=zstore_dummy( ibool(NGLLX,NGLLY,(NGLLZ+1)/2,ispec) )
+ xelm(7)=xstore_dummy( ibool(NGLLX,(NGLLY+1)/2,NGLLZ,ispec) )
+ yelm(7)=ystore_dummy( ibool(NGLLX,(NGLLY+1)/2,NGLLZ,ispec) )
+ zelm(7)=zstore_dummy( ibool(NGLLX,(NGLLY+1)/2,NGLLZ,ispec) )
+ xelm(8)=xstore_dummy( ibool(NGLLX,1,(NGLLZ+1)/2,ispec) )
+ yelm(8)=ystore_dummy( ibool(NGLLX,1,(NGLLZ+1)/2,ispec) )
+ zelm(8)=zstore_dummy( ibool(NGLLX,1,(NGLLZ+1)/2,ispec) )
+ xelm(9)=xstore_dummy( ibool(NGLLX,(NGLLY+1)/2,(NGLLZ+1)/2,ispec) )
+ yelm(9)=ystore_dummy( ibool(NGLLX,(NGLLY+1)/2,(NGLLZ+1)/2,ispec) )
+ zelm(9)=zstore_dummy( ibool(NGLLX,(NGLLY+1)/2,(NGLLZ+1)/2,ispec) )
+ endif
+
call compute_jacobian_2D_face(myrank,xelm,yelm,zelm, &
dershape2D_x,wgllwgll_yz, &
jacobian2Dw_face,normal_face,NGLLY,NGLLZ)
@@ -118,6 +154,24 @@
yelm(4)=ystore_dummy( ibool(1,1,NGLLZ,ispec) )
zelm(4)=zstore_dummy( ibool(1,1,NGLLZ,ispec) )
+ if(NGNOD2D_value == 9) then
+ xelm(5)=xstore_dummy( ibool((NGLLX+1)/2,1,1,ispec) )
+ yelm(5)=ystore_dummy( ibool((NGLLX+1)/2,1,1,ispec) )
+ zelm(5)=zstore_dummy( ibool((NGLLX+1)/2,1,1,ispec) )
+ xelm(6)=xstore_dummy( ibool(NGLLX,1,(NGLLZ+1)/2,ispec) )
+ yelm(6)=ystore_dummy( ibool(NGLLX,1,(NGLLZ+1)/2,ispec) )
+ zelm(6)=zstore_dummy( ibool(NGLLX,1,(NGLLZ+1)/2,ispec) )
+ xelm(7)=xstore_dummy( ibool((NGLLX+1)/2,1,NGLLZ,ispec) )
+ yelm(7)=ystore_dummy( ibool((NGLLX+1)/2,1,NGLLZ,ispec) )
+ zelm(7)=zstore_dummy( ibool((NGLLX+1)/2,1,NGLLZ,ispec) )
+ xelm(8)=xstore_dummy( ibool(1,1,(NGLLZ+1)/2,ispec) )
+ yelm(8)=ystore_dummy( ibool(1,1,(NGLLZ+1)/2,ispec) )
+ zelm(8)=zstore_dummy( ibool(1,1,(NGLLZ+1)/2,ispec) )
+ xelm(9)=xstore_dummy( ibool((NGLLX+1)/2,1,(NGLLZ+1)/2,ispec) )
+ yelm(9)=ystore_dummy( ibool((NGLLX+1)/2,1,(NGLLZ+1)/2,ispec) )
+ zelm(9)=zstore_dummy( ibool((NGLLX+1)/2,1,(NGLLZ+1)/2,ispec) )
+ endif
+
call compute_jacobian_2D_face(myrank,xelm,yelm,zelm, &
dershape2D_y,wgllwgll_xz, &
jacobian2Dw_face,normal_face,NGLLX,NGLLZ)
@@ -137,6 +191,24 @@
yelm(4)=ystore_dummy( ibool(1,NGLLY,NGLLZ,ispec) )
zelm(4)=zstore_dummy( ibool(1,NGLLY,NGLLZ,ispec) )
+ if(NGNOD2D_value == 9) then
+ xelm(5)=xstore_dummy( ibool((NGLLX+1)/2,NGLLY,1,ispec) )
+ yelm(5)=ystore_dummy( ibool((NGLLX+1)/2,NGLLY,1,ispec) )
+ zelm(5)=zstore_dummy( ibool((NGLLX+1)/2,NGLLY,1,ispec) )
+ xelm(6)=xstore_dummy( ibool(NGLLX,NGLLY,(NGLLZ+1)/2,ispec) )
+ yelm(6)=ystore_dummy( ibool(NGLLX,NGLLY,(NGLLZ+1)/2,ispec) )
+ zelm(6)=zstore_dummy( ibool(NGLLX,NGLLY,(NGLLZ+1)/2,ispec) )
+ xelm(7)=xstore_dummy( ibool((NGLLX+1)/2,NGLLY,NGLLZ,ispec) )
+ yelm(7)=ystore_dummy( ibool((NGLLX+1)/2,NGLLY,NGLLZ,ispec) )
+ zelm(7)=zstore_dummy( ibool((NGLLX+1)/2,NGLLY,NGLLZ,ispec) )
+ xelm(8)=xstore_dummy( ibool(1,NGLLY,(NGLLZ+1)/2,ispec) )
+ yelm(8)=ystore_dummy( ibool(1,NGLLY,(NGLLZ+1)/2,ispec) )
+ zelm(8)=zstore_dummy( ibool(1,NGLLY,(NGLLZ+1)/2,ispec) )
+ xelm(9)=xstore_dummy( ibool((NGLLX+1)/2,NGLLY,(NGLLZ+1)/2,ispec) )
+ yelm(9)=ystore_dummy( ibool((NGLLX+1)/2,NGLLY,(NGLLZ+1)/2,ispec) )
+ zelm(9)=zstore_dummy( ibool((NGLLX+1)/2,NGLLY,(NGLLZ+1)/2,ispec) )
+ endif
+
call compute_jacobian_2D_face(myrank,xelm,yelm,zelm, &
dershape2D_y, wgllwgll_xz, &
jacobian2Dw_face,normal_face,NGLLX,NGLLZ)
@@ -157,6 +229,24 @@
yelm(4)=ystore_dummy( ibool(1,NGLLY,1,ispec) )
zelm(4)=zstore_dummy( ibool(1,NGLLY,1,ispec) )
+ if(NGNOD2D_value == 9) then
+ xelm(5)=xstore_dummy( ibool((NGLLX+1)/2,1,1,ispec) )
+ yelm(5)=ystore_dummy( ibool((NGLLX+1)/2,1,1,ispec) )
+ zelm(5)=zstore_dummy( ibool((NGLLX+1)/2,1,1,ispec) )
+ xelm(6)=xstore_dummy( ibool(NGLLX,(NGLLY+1)/2,1,ispec) )
+ yelm(6)=ystore_dummy( ibool(NGLLX,(NGLLY+1)/2,1,ispec) )
+ zelm(6)=zstore_dummy( ibool(NGLLX,(NGLLY+1)/2,1,ispec) )
+ xelm(7)=xstore_dummy( ibool((NGLLX+1)/2,NGLLY,1,ispec) )
+ yelm(7)=ystore_dummy( ibool((NGLLX+1)/2,NGLLY,1,ispec) )
+ zelm(7)=zstore_dummy( ibool((NGLLX+1)/2,NGLLY,1,ispec) )
+ xelm(8)=xstore_dummy( ibool(1,(NGLLY+1)/2,1,ispec) )
+ yelm(8)=ystore_dummy( ibool(1,(NGLLY+1)/2,1,ispec) )
+ zelm(8)=zstore_dummy( ibool(1,(NGLLY+1)/2,1,ispec) )
+ xelm(9)=xstore_dummy( ibool((NGLLX+1)/2,(NGLLY+1)/2,1,ispec) )
+ yelm(9)=ystore_dummy( ibool((NGLLX+1)/2,(NGLLY+1)/2,1,ispec) )
+ zelm(9)=zstore_dummy( ibool((NGLLX+1)/2,(NGLLY+1)/2,1,ispec) )
+ endif
+
call compute_jacobian_2D_face(myrank,xelm,yelm,zelm,&
dershape2D_bottom,wgllwgll_xy, &
jacobian2Dw_face,normal_face,NGLLX,NGLLY)
@@ -176,6 +266,24 @@
yelm(4)=ystore_dummy( ibool(1,NGLLY,NGLLZ,ispec) )
zelm(4)=zstore_dummy( ibool(1,NGLLY,NGLLZ,ispec) )
+ if(NGNOD2D_value == 9) then
+ xelm(5)=xstore_dummy( ibool((NGLLX+1)/2,1,NGLLZ,ispec) )
+ yelm(5)=ystore_dummy( ibool((NGLLX+1)/2,1,NGLLZ,ispec) )
+ zelm(5)=zstore_dummy( ibool((NGLLX+1)/2,1,NGLLZ,ispec) )
+ xelm(6)=xstore_dummy( ibool(NGLLX,(NGLLY+1)/2,NGLLZ,ispec) )
+ yelm(6)=ystore_dummy( ibool(NGLLX,(NGLLY+1)/2,NGLLZ,ispec) )
+ zelm(6)=zstore_dummy( ibool(NGLLX,(NGLLY+1)/2,NGLLZ,ispec) )
+ xelm(7)=xstore_dummy( ibool((NGLLX+1)/2,NGLLY,NGLLZ,ispec) )
+ yelm(7)=ystore_dummy( ibool((NGLLX+1)/2,NGLLY,NGLLZ,ispec) )
+ zelm(7)=zstore_dummy( ibool((NGLLX+1)/2,NGLLY,NGLLZ,ispec) )
+ xelm(8)=xstore_dummy( ibool(1,(NGLLY+1)/2,NGLLZ,ispec) )
+ yelm(8)=ystore_dummy( ibool(1,(NGLLY+1)/2,NGLLZ,ispec) )
+ zelm(8)=zstore_dummy( ibool(1,(NGLLY+1)/2,NGLLZ,ispec) )
+ xelm(9)=xstore_dummy( ibool((NGLLX+1)/2,(NGLLY+1)/2,NGLLZ,ispec) )
+ yelm(9)=ystore_dummy( ibool((NGLLX+1)/2,(NGLLY+1)/2,NGLLZ,ispec) )
+ zelm(9)=zstore_dummy( ibool((NGLLX+1)/2,(NGLLY+1)/2,NGLLZ,ispec) )
+ endif
+
call compute_jacobian_2D_face(myrank,xelm,yelm,zelm,&
dershape2D_top, wgllwgll_xy, &
jacobian2Dw_face,normal_face,NGLLX,NGLLY)
@@ -243,7 +351,7 @@
! distinguish if single or double precision for reals
if(CUSTOM_REAL == SIZE_REAL) then
- jacobian2Dw_face(i,j) = sngl(jacobian * wgllwgll(i,j) )
+ jacobian2Dw_face(i,j) = sngl(jacobian * wgllwgll(i,j))
normal_face(1,i,j)=sngl(unx/jacobian)
normal_face(2,i,j)=sngl(uny/jacobian)
normal_face(3,i,j)=sngl(unz/jacobian)
@@ -254,37 +362,6 @@
normal_face(3,i,j)=unz/jacobian
endif
- !debug
- ! note: normal values could be almost zero and lead to underflow errors
- !if(FIX_UNDERFLOW_PROBLEM) then
- ! if( abs(normal_face(1,i,j)) < VERYSMALLVAL ) &
- ! normal_face(1,i,j) = sign(1.0_CUSTOM_REAL,normal_face(1,i,j)) * VERYSMALLVAL
- ! if( abs(normal_face(2,i,j)) < VERYSMALLVAL ) &
- ! normal_face(2,i,j) = sign(1.0_CUSTOM_REAL,normal_face(2,i,j)) * VERYSMALLVAL
- ! if( abs(normal_face(3,i,j)) < VERYSMALLVAL ) &
- ! normal_face(3,i,j) = sign(1.0_CUSTOM_REAL,normal_face(3,i,j)) * VERYSMALLVAL
- !endif
- ! re-calculates length of normal (might have rounding errors)
- !jacobian = sqrt(dble(normal_face(1,i,j))**2 + dble(normal_face(2,i,j))**2 + dble(normal_face(3,i,j))**2)
- !if( abs(jacobian - 1.d0) > TINYVAL ) then
- ! if(jacobian <= ZERO) call exit_MPI(myrank,'normal length undefined')
- ! ! re-normalizes
- ! normal_face(:,i,j) = normal_face(:,i,j) / jacobian
- !endif
- !debug
- !if( isNan(normal_face(1,i,j)) ) then
- ! print*,'error normal_face 1:',normal_face(:,i,j)
- ! stop 'error normal_face'
- !endif
- !if( isNan(normal_face(2,i,j)) ) then
- ! print*,'error normal_face 2:',normal_face(:,i,j)
- ! stop 'error normal_face'
- !endif
- !if( isNan(normal_face(3,i,j)) ) then
- ! print*,'error normal_face 3:',normal_face(:,i,j)
- ! stop 'error normal_face'
- !endif
-
enddo
enddo
Modified: seismo/3D/SPECFEM3D/trunk/src/shared/get_shape2D.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/shared/get_shape2D.f90 2012-10-13 18:42:30 UTC (rev 20836)
+++ seismo/3D/SPECFEM3D/trunk/src/shared/get_shape2D.f90 2012-10-13 22:18:10 UTC (rev 20837)
@@ -234,4 +234,3 @@
end subroutine get_shape2D_9
-
Modified: seismo/3D/SPECFEM3D/trunk/src/shared/get_shape3D.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/shared/get_shape3D.f90 2012-10-13 18:42:30 UTC (rev 20836)
+++ seismo/3D/SPECFEM3D/trunk/src/shared/get_shape3D.f90 2012-10-13 22:18:10 UTC (rev 20837)
@@ -55,14 +55,12 @@
double precision, parameter :: ONE_EIGHTH = 0.125d0
! check that the parameter file is correct
- if(NGNOD /= 8) call exit_MPI(myrank,'elements should have 8 control nodes')
+ if(NGNOD /= 8 .and. NGNOD /= 27) call exit_MPI(myrank,'elements should have 8 or 27 control nodes')
! ***
! *** create 3D shape functions and jacobian
! ***
-!--- case of a 3D 8-node element (Dhatt-Touzot p. 115)
-
do i=1,NGLLX
do j=1,NGLLY
do k=1,NGLLZ
@@ -71,6 +69,9 @@
eta = yigll(j)
gamma = zigll(k)
+!--- case of a 3D 8-node element (Dhatt-Touzot p. 115)
+ if(NGNOD == 8) then
+
ra1 = one + xi
ra2 = one - xi
@@ -116,6 +117,14 @@
dershape3D(3,7,i,j,k) = ONE_EIGHTH*ra1*rb1
dershape3D(3,8,i,j,k) = ONE_EIGHTH*ra2*rb1
+ else
+
+ ! note: put further initialization for ngnod == 27 into subroutine
+ ! to avoid compilation errors in case ngnod == 8
+ call get_shape3D_27(NGNOD,shape3D,dershape3D,xi,eta,gamma,i,j,k)
+
+ endif
+
enddo
enddo
enddo
@@ -140,10 +149,10 @@
! sum of shape functions should be one
! sum of derivative of shape functions should be zero
- if(abs(sumshape-one) > TINYVAL) call exit_MPI(myrank,'error shape functions')
- if(abs(sumdershapexi) > TINYVAL) call exit_MPI(myrank,'error derivative xi shape functions')
- if(abs(sumdershapeeta) > TINYVAL) call exit_MPI(myrank,'error derivative eta shape functions')
- if(abs(sumdershapegamma) > TINYVAL) call exit_MPI(myrank,'error derivative gamma shape functions')
+ if(abs(sumshape-one) > TINYVAL) call exit_MPI(myrank,'error in 3D shape functions')
+ if(abs(sumdershapexi) > TINYVAL) call exit_MPI(myrank,'error in xi derivative of 3D shape functions')
+ if(abs(sumdershapeeta) > TINYVAL) call exit_MPI(myrank,'error in eta derivative of 3D shape functions')
+ if(abs(sumdershapegamma) > TINYVAL) call exit_MPI(myrank,'error in gamma derivative of 3D shape functions')
enddo
enddo
@@ -157,52 +166,107 @@
! 3D shape functions for given, single xi/eta/gamma location
- subroutine eval_shape3D_single(myrank,shape3D,xi,eta,gamma)
+ subroutine eval_shape3D_single(myrank,shape3D,xi,eta,gamma,NGNOD_value)
implicit none
include "constants.h"
- integer :: myrank
+ integer :: myrank,NGNOD_value
! 3D shape functions
- double precision :: shape3D(NGNOD)
+ double precision :: shape3D(NGNOD_value)
! location
double precision :: xi,eta,gamma
! local parameters
+ double precision, parameter :: ONE_EIGHTH = 0.125d0
double precision :: ra1,ra2,rb1,rb2,rc1,rc2
- double precision, parameter :: ONE_EIGHTH = 0.125d0
+ double precision :: l1xi,l2xi,l3xi,l1eta,l2eta,l3eta,l1gamma,l2gamma,l3gamma
double precision :: sumshape
integer :: ia
! check that the parameter file is correct
- if(NGNOD /= 8) call exit_MPI(myrank,'elements should have 8 control nodes')
+ if(NGNOD_value /= 8 .and. NGNOD_value /= 27) call exit_MPI(myrank,'elements should have 8 or 27 control nodes')
+ ! shape functions
+
!--- case of a 3D 8-node element (Dhatt-Touzot p. 115)
- ra1 = one + xi
- ra2 = one - xi
+ if(NGNOD_value == 8) then
- rb1 = one + eta
- rb2 = one - eta
+ ra1 = one + xi
+ ra2 = one - xi
- rc1 = one + gamma
- rc2 = one - gamma
+ rb1 = one + eta
+ rb2 = one - eta
- ! shape functions
- shape3D(1) = ONE_EIGHTH*ra2*rb2*rc2
- shape3D(2) = ONE_EIGHTH*ra1*rb2*rc2
- shape3D(3) = ONE_EIGHTH*ra1*rb1*rc2
- shape3D(4) = ONE_EIGHTH*ra2*rb1*rc2
- shape3D(5) = ONE_EIGHTH*ra2*rb2*rc1
- shape3D(6) = ONE_EIGHTH*ra1*rb2*rc1
- shape3D(7) = ONE_EIGHTH*ra1*rb1*rc1
- shape3D(8) = ONE_EIGHTH*ra2*rb1*rc1
+ rc1 = one + gamma
+ rc2 = one - gamma
+ shape3D(1) = ONE_EIGHTH*ra2*rb2*rc2
+ shape3D(2) = ONE_EIGHTH*ra1*rb2*rc2
+ shape3D(3) = ONE_EIGHTH*ra1*rb1*rc2
+ shape3D(4) = ONE_EIGHTH*ra2*rb1*rc2
+ shape3D(5) = ONE_EIGHTH*ra2*rb2*rc1
+ shape3D(6) = ONE_EIGHTH*ra1*rb2*rc1
+ shape3D(7) = ONE_EIGHTH*ra1*rb1*rc1
+ shape3D(8) = ONE_EIGHTH*ra2*rb1*rc1
+
+ else
+
+ l1xi=HALF*xi*(xi-ONE)
+ l2xi=ONE-xi**2
+ l3xi=HALF*xi*(xi+ONE)
+
+ l1eta=HALF*eta*(eta-ONE)
+ l2eta=ONE-eta**2
+ l3eta=HALF*eta*(eta+ONE)
+
+ l1gamma=HALF*gamma*(gamma-ONE)
+ l2gamma=ONE-gamma**2
+ l3gamma=HALF*gamma*(gamma+ONE)
+
+! corner nodes
+ shape3D(1)=l1xi*l1eta*l1gamma
+ shape3D(2)=l3xi*l1eta*l1gamma
+ shape3D(3)=l3xi*l3eta*l1gamma
+ shape3D(4)=l1xi*l3eta*l1gamma
+ shape3D(5)=l1xi*l1eta*l3gamma
+ shape3D(6)=l3xi*l1eta*l3gamma
+ shape3D(7)=l3xi*l3eta*l3gamma
+ shape3D(8)=l1xi*l3eta*l3gamma
+
+! midside nodes
+ shape3D(9)=l2xi*l1eta*l1gamma
+ shape3D(10)=l3xi*l2eta*l1gamma
+ shape3D(11)=l2xi*l3eta*l1gamma
+ shape3D(12)=l1xi*l2eta*l1gamma
+ shape3D(13)=l1xi*l1eta*l2gamma
+ shape3D(14)=l3xi*l1eta*l2gamma
+ shape3D(15)=l3xi*l3eta*l2gamma
+ shape3D(16)=l1xi*l3eta*l2gamma
+ shape3D(17)=l2xi*l1eta*l3gamma
+ shape3D(18)=l3xi*l2eta*l3gamma
+ shape3D(19)=l2xi*l3eta*l3gamma
+ shape3D(20)=l1xi*l2eta*l3gamma
+
+! side center nodes
+ shape3D(21)=l2xi*l2eta*l1gamma
+ shape3D(22)=l2xi*l1eta*l2gamma
+ shape3D(23)=l3xi*l2eta*l2gamma
+ shape3D(24)=l2xi*l3eta*l2gamma
+ shape3D(25)=l1xi*l2eta*l2gamma
+ shape3D(26)=l2xi*l2eta*l3gamma
+
+! center node
+ shape3D(27)=l2xi*l2eta*l2gamma
+
+ endif
+
! check the shape functions
sumshape = ZERO
- do ia=1,NGNOD
+ do ia=1,NGNOD_value
sumshape = sumshape + shape3D(ia)
enddo
@@ -267,3 +331,178 @@
end subroutine eval_shape3D_element_corners
+!
+!-------------------------------------------------------------------------------------------------
+!
+
+!--- case of a 3D 27-node element
+
+ subroutine get_shape3D_27(NGNOD_value,shape3D,dershape3D,xi,eta,gamma,i,j,k)
+
+ implicit none
+
+ include "constants.h"
+
+ integer :: NGNOD_value,i,j,k
+
+! 3D shape functions and their derivatives
+ double precision shape3D(NGNOD_value,NGLLX,NGLLY,NGLLZ)
+ double precision dershape3D(NDIM,NGNOD_value,NGLLX,NGLLY,NGLLZ)
+
+! location of the nodes of the 3D quadrilateral elements
+ double precision xi,eta,gamma
+ double precision l1xi,l2xi,l3xi,l1eta,l2eta,l3eta,l1gamma,l2gamma,l3gamma
+ double precision l1pxi,l2pxi,l3pxi,l1peta,l2peta,l3peta,l1pgamma,l2pgamma,l3pgamma
+
+ l1xi=HALF*xi*(xi-ONE)
+ l2xi=ONE-xi**2
+ l3xi=HALF*xi*(xi+ONE)
+
+ l1pxi=xi-HALF
+ l2pxi=-TWO*xi
+ l3pxi=xi+HALF
+
+ l1eta=HALF*eta*(eta-ONE)
+ l2eta=ONE-eta**2
+ l3eta=HALF*eta*(eta+ONE)
+
+ l1peta=eta-HALF
+ l2peta=-TWO*eta
+ l3peta=eta+HALF
+
+ l1gamma=HALF*gamma*(gamma-ONE)
+ l2gamma=ONE-gamma**2
+ l3gamma=HALF*gamma*(gamma+ONE)
+
+ l1pgamma=gamma-HALF
+ l2pgamma=-TWO*gamma
+ l3pgamma=gamma+HALF
+
+! corner nodes
+ shape3D(1,i,j,k)=l1xi*l1eta*l1gamma
+ shape3D(2,i,j,k)=l3xi*l1eta*l1gamma
+ shape3D(3,i,j,k)=l3xi*l3eta*l1gamma
+ shape3D(4,i,j,k)=l1xi*l3eta*l1gamma
+ shape3D(5,i,j,k)=l1xi*l1eta*l3gamma
+ shape3D(6,i,j,k)=l3xi*l1eta*l3gamma
+ shape3D(7,i,j,k)=l3xi*l3eta*l3gamma
+ shape3D(8,i,j,k)=l1xi*l3eta*l3gamma
+
+ dershape3D(1,1,i,j,k)=l1pxi*l1eta*l1gamma
+ dershape3D(1,2,i,j,k)=l3pxi*l1eta*l1gamma
+ dershape3D(1,3,i,j,k)=l3pxi*l3eta*l1gamma
+ dershape3D(1,4,i,j,k)=l1pxi*l3eta*l1gamma
+ dershape3D(1,5,i,j,k)=l1pxi*l1eta*l3gamma
+ dershape3D(1,6,i,j,k)=l3pxi*l1eta*l3gamma
+ dershape3D(1,7,i,j,k)=l3pxi*l3eta*l3gamma
+ dershape3D(1,8,i,j,k)=l1pxi*l3eta*l3gamma
+
+ dershape3D(2,1,i,j,k)=l1xi*l1peta*l1gamma
+ dershape3D(2,2,i,j,k)=l3xi*l1peta*l1gamma
+ dershape3D(2,3,i,j,k)=l3xi*l3peta*l1gamma
+ dershape3D(2,4,i,j,k)=l1xi*l3peta*l1gamma
+ dershape3D(2,5,i,j,k)=l1xi*l1peta*l3gamma
+ dershape3D(2,6,i,j,k)=l3xi*l1peta*l3gamma
+ dershape3D(2,7,i,j,k)=l3xi*l3peta*l3gamma
+ dershape3D(2,8,i,j,k)=l1xi*l3peta*l3gamma
+
+ dershape3D(3,1,i,j,k)=l1xi*l1eta*l1pgamma
+ dershape3D(3,2,i,j,k)=l3xi*l1eta*l1pgamma
+ dershape3D(3,3,i,j,k)=l3xi*l3eta*l1pgamma
+ dershape3D(3,4,i,j,k)=l1xi*l3eta*l1pgamma
+ dershape3D(3,5,i,j,k)=l1xi*l1eta*l3pgamma
+ dershape3D(3,6,i,j,k)=l3xi*l1eta*l3pgamma
+ dershape3D(3,7,i,j,k)=l3xi*l3eta*l3pgamma
+ dershape3D(3,8,i,j,k)=l1xi*l3eta*l3pgamma
+
+! midside nodes
+ shape3D(9,i,j,k)=l2xi*l1eta*l1gamma
+ shape3D(10,i,j,k)=l3xi*l2eta*l1gamma
+ shape3D(11,i,j,k)=l2xi*l3eta*l1gamma
+ shape3D(12,i,j,k)=l1xi*l2eta*l1gamma
+ shape3D(13,i,j,k)=l1xi*l1eta*l2gamma
+ shape3D(14,i,j,k)=l3xi*l1eta*l2gamma
+ shape3D(15,i,j,k)=l3xi*l3eta*l2gamma
+ shape3D(16,i,j,k)=l1xi*l3eta*l2gamma
+ shape3D(17,i,j,k)=l2xi*l1eta*l3gamma
+ shape3D(18,i,j,k)=l3xi*l2eta*l3gamma
+ shape3D(19,i,j,k)=l2xi*l3eta*l3gamma
+ shape3D(20,i,j,k)=l1xi*l2eta*l3gamma
+
+ dershape3D(1,9,i,j,k)=l2pxi*l1eta*l1gamma
+ dershape3D(1,10,i,j,k)=l3pxi*l2eta*l1gamma
+ dershape3D(1,11,i,j,k)=l2pxi*l3eta*l1gamma
+ dershape3D(1,12,i,j,k)=l1pxi*l2eta*l1gamma
+ dershape3D(1,13,i,j,k)=l1pxi*l1eta*l2gamma
+ dershape3D(1,14,i,j,k)=l3pxi*l1eta*l2gamma
+ dershape3D(1,15,i,j,k)=l3pxi*l3eta*l2gamma
+ dershape3D(1,16,i,j,k)=l1pxi*l3eta*l2gamma
+ dershape3D(1,17,i,j,k)=l2pxi*l1eta*l3gamma
+ dershape3D(1,18,i,j,k)=l3pxi*l2eta*l3gamma
+ dershape3D(1,19,i,j,k)=l2pxi*l3eta*l3gamma
+ dershape3D(1,20,i,j,k)=l1pxi*l2eta*l3gamma
+
+ dershape3D(2,9,i,j,k)=l2xi*l1peta*l1gamma
+ dershape3D(2,10,i,j,k)=l3xi*l2peta*l1gamma
+ dershape3D(2,11,i,j,k)=l2xi*l3peta*l1gamma
+ dershape3D(2,12,i,j,k)=l1xi*l2peta*l1gamma
+ dershape3D(2,13,i,j,k)=l1xi*l1peta*l2gamma
+ dershape3D(2,14,i,j,k)=l3xi*l1peta*l2gamma
+ dershape3D(2,15,i,j,k)=l3xi*l3peta*l2gamma
+ dershape3D(2,16,i,j,k)=l1xi*l3peta*l2gamma
+ dershape3D(2,17,i,j,k)=l2xi*l1peta*l3gamma
+ dershape3D(2,18,i,j,k)=l3xi*l2peta*l3gamma
+ dershape3D(2,19,i,j,k)=l2xi*l3peta*l3gamma
+ dershape3D(2,20,i,j,k)=l1xi*l2peta*l3gamma
+
+ dershape3D(3,9,i,j,k)=l2xi*l1eta*l1pgamma
+ dershape3D(3,10,i,j,k)=l3xi*l2eta*l1pgamma
+ dershape3D(3,11,i,j,k)=l2xi*l3eta*l1pgamma
+ dershape3D(3,12,i,j,k)=l1xi*l2eta*l1pgamma
+ dershape3D(3,13,i,j,k)=l1xi*l1eta*l2pgamma
+ dershape3D(3,14,i,j,k)=l3xi*l1eta*l2pgamma
+ dershape3D(3,15,i,j,k)=l3xi*l3eta*l2pgamma
+ dershape3D(3,16,i,j,k)=l1xi*l3eta*l2pgamma
+ dershape3D(3,17,i,j,k)=l2xi*l1eta*l3pgamma
+ dershape3D(3,18,i,j,k)=l3xi*l2eta*l3pgamma
+ dershape3D(3,19,i,j,k)=l2xi*l3eta*l3pgamma
+ dershape3D(3,20,i,j,k)=l1xi*l2eta*l3pgamma
+
+! side center nodes
+ shape3D(21,i,j,k)=l2xi*l2eta*l1gamma
+ shape3D(22,i,j,k)=l2xi*l1eta*l2gamma
+ shape3D(23,i,j,k)=l3xi*l2eta*l2gamma
+ shape3D(24,i,j,k)=l2xi*l3eta*l2gamma
+ shape3D(25,i,j,k)=l1xi*l2eta*l2gamma
+ shape3D(26,i,j,k)=l2xi*l2eta*l3gamma
+
+ dershape3D(1,21,i,j,k)=l2pxi*l2eta*l1gamma
+ dershape3D(1,22,i,j,k)=l2pxi*l1eta*l2gamma
+ dershape3D(1,23,i,j,k)=l3pxi*l2eta*l2gamma
+ dershape3D(1,24,i,j,k)=l2pxi*l3eta*l2gamma
+ dershape3D(1,25,i,j,k)=l1pxi*l2eta*l2gamma
+ dershape3D(1,26,i,j,k)=l2pxi*l2eta*l3gamma
+
+ dershape3D(2,21,i,j,k)=l2xi*l2peta*l1gamma
+ dershape3D(2,22,i,j,k)=l2xi*l1peta*l2gamma
+ dershape3D(2,23,i,j,k)=l3xi*l2peta*l2gamma
+ dershape3D(2,24,i,j,k)=l2xi*l3peta*l2gamma
+ dershape3D(2,25,i,j,k)=l1xi*l2peta*l2gamma
+ dershape3D(2,26,i,j,k)=l2xi*l2peta*l3gamma
+
+ dershape3D(3,21,i,j,k)=l2xi*l2eta*l1pgamma
+ dershape3D(3,22,i,j,k)=l2xi*l1eta*l2pgamma
+ dershape3D(3,23,i,j,k)=l3xi*l2eta*l2pgamma
+ dershape3D(3,24,i,j,k)=l2xi*l3eta*l2pgamma
+ dershape3D(3,25,i,j,k)=l1xi*l2eta*l2pgamma
+ dershape3D(3,26,i,j,k)=l2xi*l2eta*l3pgamma
+
+! center node
+ shape3D(27,i,j,k)=l2xi*l2eta*l2gamma
+
+ dershape3D(1,27,i,j,k)=l2pxi*l2eta*l2gamma
+ dershape3D(2,27,i,j,k)=l2xi*l2peta*l2gamma
+ dershape3D(3,27,i,j,k)=l2xi*l2eta*l2pgamma
+
+ end subroutine get_shape3D_27
+
Modified: seismo/3D/SPECFEM3D/trunk/src/shared/prepare_assemble_MPI.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/shared/prepare_assemble_MPI.f90 2012-10-13 18:42:30 UTC (rev 20836)
+++ seismo/3D/SPECFEM3D/trunk/src/shared/prepare_assemble_MPI.f90 2012-10-13 22:18:10 UTC (rev 20837)
@@ -42,7 +42,7 @@
! spectral element indexing
! ( nelmnts = number of spectral elements
-! NGNOD = number of element corners (8)
+! NGNOD = number of element control points
! knods = corner indices array )
integer, intent(in) :: nelmnts
integer, dimension(NGNOD,nelmnts), intent(in) :: knods
@@ -69,7 +69,7 @@
logical, dimension(:),allocatable :: mask_ibool_asteroid
integer :: ixmin, ixmax, iymin, iymax, izmin, izmax
- integer, dimension(NGNOD) :: n
+ integer, dimension(NGNOD_EIGHT_CORNERS) :: n
integer :: e1, e2, e3, e4
integer :: ispec,k,ix,iy,iz,ier,itype,iglob
integer :: npoin_interface_asteroid
@@ -92,7 +92,7 @@
! type of interface: (1) corner point, (2) edge, (4) face
itype = my_interfaces(2,ispec_interface,num_interface)
! gets spectral element corner indices (defines all nodes of face/edge)
- do k = 1, NGNOD
+ do k = 1, NGNOD_EIGHT_CORNERS
n(k) = knods(k,ispec)
end do
@@ -151,7 +151,7 @@
include "constants.h"
! corner node indices per spectral element (8)
- integer, dimension(NGNOD), intent(in) :: n
+ integer, dimension(NGNOD_EIGHT_CORNERS), intent(in) :: n
! interface type & nodes
integer, intent(in) :: itype, e1, e2, e3, e4
Modified: seismo/3D/SPECFEM3D/trunk/src/shared/read_parameter_file.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/shared/read_parameter_file.f90 2012-10-13 18:42:30 UTC (rev 20836)
+++ seismo/3D/SPECFEM3D/trunk/src/shared/read_parameter_file.f90 2012-10-13 22:18:10 UTC (rev 20837)
@@ -248,7 +248,7 @@
! one cannot use a Heaviside source for the movies
if((MOVIE_SURFACE .or. MOVIE_VOLUME) .and. sqrt(minval_hdur**2 + HDUR_MOVIE**2) < TINYVAL) &
stop 'hdur too small for movie creation, movies do not make sense for Heaviside source'
- endif
+ endif
! converts all string characters to lowercase
irange = iachar('a') - iachar('A')
Modified: seismo/3D/SPECFEM3D/trunk/src/shared/save_header_file.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/shared/save_header_file.f90 2012-10-13 18:42:30 UTC (rev 20836)
+++ seismo/3D/SPECFEM3D/trunk/src/shared/save_header_file.f90 2012-10-13 22:18:10 UTC (rev 20837)
@@ -44,14 +44,14 @@
character(len=256) HEADER_FILE
integer nfaces_surface_glob_ext_mesh
-
+
NAMELIST/MESHER/ABSORB_FREE_SURFACE_VAL
if (ABSORB_INSTEAD_OF_FREE_SURFACE) then
ABSORB_FREE_SURFACE_VAL = .true.
else
ABSORB_FREE_SURFACE_VAL = .false.
- endif
+ endif
! copy number of elements and points in an include file for the solver
call get_value_string(HEADER_FILE, 'solver.HEADER_FILE', &
@@ -125,7 +125,7 @@
write(IOUT,*) '! (if significantly less, you waste a significant amount of memory)'
write(IOUT,*) '!'
write(IOUT,*) '! check parameter to ensure the code has been compiled with the right values:'
- write(IOUT,NML=MESHER)
+ write(IOUT,NML=MESHER)
write(IOUT,*)
close(IOUT)
Modified: seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_add_sources_elastic.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_add_sources_elastic.f90 2012-10-13 18:42:30 UTC (rev 20836)
+++ seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_add_sources_elastic.f90 2012-10-13 22:18:10 UTC (rev 20837)
@@ -142,7 +142,7 @@
! change dt -> DT
call compute_add_sources_el_cuda(Mesh_pointer, phase_is_inner, &
NSOURCES, stf_pre_compute, myrank)
-
+
endif
else ! .NOT. GPU_MODE
@@ -177,7 +177,7 @@
! This is the expression of a Ricker; should be changed according maybe to the Par_file.
stf = comp_source_time_function_rickr(dble(it-1)*DT-t0-tshift_cmt(isource),f0)
- ! add the inclined force source array
+ ! add the inclined force source array
! distinguish between single and double precision for reals
if(CUSTOM_REAL == SIZE_REAL) then
stf_used = sngl(stf)
@@ -193,7 +193,7 @@
enddo
enddo
enddo
-
+
else
if( USE_RICKER_IPATI) then
@@ -463,7 +463,7 @@
! This is the expression of a Ricker; should be changed according maybe to the Par_file.
stf = comp_source_time_function_rickr(dble(NSTEP-it)*DT-t0-tshift_cmt(isource),f0)
- ! add the inclined force source array
+ ! add the inclined force source array
! distinguish between single and double precision for reals
if(CUSTOM_REAL == SIZE_REAL) then
stf_used = sngl(stf)
Modified: seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_add_sources_poroelastic.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_add_sources_poroelastic.f90 2012-10-13 18:42:30 UTC (rev 20836)
+++ seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_add_sources_poroelastic.f90 2012-10-13 22:18:10 UTC (rev 20837)
@@ -157,7 +157,7 @@
stf = comp_source_time_function_rickr(dble(it-1)*DT-t0-tshift_cmt(isource),f0)
!stf_used = comp_source_time_function_gauss(dble(it-1)*DT-t0-tshift_cmt(isource),f0)
- ! add the inclined force source array
+ ! add the inclined force source array
! the source is applied to both solid and fluid phase: bulk source.
! distinguish between single and double precision for reals
@@ -176,7 +176,7 @@
(1._CUSTOM_REAL - phil/tortl) * sourcearrays(isource,:,i,j,k) * stf_used
! fluid phase
accelw(:,iglob) = accelw(:,iglob) + &
- (1._CUSTOM_REAL - rhol_f/rhol_bar) * sourcearrays(isource,:,i,j,k) * stf_used
+ (1._CUSTOM_REAL - rhol_f/rhol_bar) * sourcearrays(isource,:,i,j,k) * stf_used
enddo
enddo
enddo
@@ -417,7 +417,7 @@
! should be the same than for the forward simulation (check above)
stf = comp_source_time_function_rickr(dble(NSTEP-it)*DT-t0-tshift_cmt(isource),f0)
- ! add the inclined force source array
+ ! add the inclined force source array
! the source is applied to both solid and fluid phase: bulk source
! note: time step is now at NSTEP-it
Modified: seismo/3D/SPECFEM3D/trunk/src/specfem3D/initialize_simulation.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/specfem3D/initialize_simulation.f90 2012-10-13 18:42:30 UTC (rev 20836)
+++ seismo/3D/SPECFEM3D/trunk/src/specfem3D/initialize_simulation.f90 2012-10-13 22:18:10 UTC (rev 20837)
@@ -281,7 +281,7 @@
OUTPUT_FILES_PATH(1:len_trim(OUTPUT_FILES_PATH))//'/values_from_mesher.h')
open(unit=IOUT,file=HEADER_FILE,status='old')
- read(IOUT,NML=MESHER)
+ read(IOUT,NML=MESHER)
close(IOUT)
if (ABSORB_INSTEAD_OF_FREE_SURFACE .NEQV. ABSORB_FREE_SURFACE_VAL) then
Modified: seismo/3D/SPECFEM3D/trunk/src/specfem3D/locate_source.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/specfem3D/locate_source.f90 2012-10-13 18:42:30 UTC (rev 20836)
+++ seismo/3D/SPECFEM3D/trunk/src/specfem3D/locate_source.f90 2012-10-13 22:18:10 UTC (rev 20837)
@@ -135,7 +135,7 @@
double precision, dimension(NSOURCES) :: lat,long,depth
- double precision, dimension(NSOURCES) :: factor_force_source
+ double precision, dimension(NSOURCES) :: factor_force_source
double precision, dimension(NSOURCES) :: comp_dir_vect_source_E
double precision, dimension(NSOURCES) :: comp_dir_vect_source_N
double precision, dimension(NSOURCES) :: comp_dir_vect_source_Z_UP
Modified: seismo/3D/SPECFEM3D/trunk/src/specfem3D/model_update.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/specfem3D/model_update.f90 2012-10-13 18:42:30 UTC (rev 20836)
+++ seismo/3D/SPECFEM3D/trunk/src/specfem3D/model_update.f90 2012-10-13 22:18:10 UTC (rev 20837)
@@ -867,7 +867,7 @@
rmass_acoustic_new = 0._CUSTOM_REAL
rmass_solid_poroelastic_new = 0._CUSTOM_REAL
rmass_fluid_poroelastic_new = 0._CUSTOM_REAL
-
+
!-------- attenuation -------
! read the proc*attenuation.vtk for the old model in LOCAL_PATH and store qmu_attenuation_store
Modified: seismo/3D/SPECFEM3D/trunk/src/specfem3D/noise_tomography.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/specfem3D/noise_tomography.f90 2012-10-13 18:42:30 UTC (rev 20836)
+++ seismo/3D/SPECFEM3D/trunk/src/specfem3D/noise_tomography.f90 2012-10-13 22:18:10 UTC (rev 20837)
@@ -455,7 +455,7 @@
if(NSPEC_TOP .gt. 0) then
! check integer size limit: size of b_reclen_field must fit onto an 4-byte integer
- if( NSPEC_TOP > 2147483647 / (CUSTOM_REAL * NGLLSQUARE * NDIM) ) then
+ if( NSPEC_TOP > 2147483646 / (CUSTOM_REAL * NGLLSQUARE * NDIM) ) then
print *,'reclen of noise surface_movie needed exceeds integer 4-byte limit: ',reclen
print *,' ',CUSTOM_REAL, NDIM, NGLLSQUARE, NSPEC_TOP
print*,'bit size fortran: ',bit_size(NSPEC_TOP)
Modified: seismo/3D/SPECFEM3D/trunk/src/specfem3D/prepare_timerun.F90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/specfem3D/prepare_timerun.F90 2012-10-13 18:42:30 UTC (rev 20836)
+++ seismo/3D/SPECFEM3D/trunk/src/specfem3D/prepare_timerun.F90 2012-10-13 22:18:10 UTC (rev 20837)
@@ -107,9 +107,9 @@
! total time per time step:
! T_total = dt * nspec
!
- ! total time using nproc processes (slices) for NSTEP time steps:
+ ! total time using nproc processes (slices) for NSTEP time steps:
! T_simulation = T_total * NSTEP / nproc
-
+
endif
! synchronize all the processes
@@ -801,7 +801,7 @@
b_reclen_field = CUSTOM_REAL * NDIM * NGLLSQUARE * num_abs_boundary_faces
! check integer size limit: size of b_reclen_field must fit onto an 4-byte integer
- if( num_abs_boundary_faces > 2147483647 / (CUSTOM_REAL * NDIM * NGLLSQUARE) ) then
+ if( num_abs_boundary_faces > 2147483646 / (CUSTOM_REAL * NDIM * NGLLSQUARE) ) then
print *,'reclen needed exceeds integer 4-byte limit: ',b_reclen_field
print *,' ',CUSTOM_REAL, NDIM, NGLLSQUARE, num_abs_boundary_faces
print*,'bit size fortran: ',bit_size(b_reclen_field)
@@ -850,7 +850,7 @@
b_reclen_potential = CUSTOM_REAL * NGLLSQUARE * num_abs_boundary_faces
! check integer size limit: size of b_reclen_potential must fit onto an 4-byte integer
- if( num_abs_boundary_faces > 2147483647 / (CUSTOM_REAL * NGLLSQUARE) ) then
+ if( num_abs_boundary_faces > 2147483646 / (CUSTOM_REAL * NGLLSQUARE) ) then
print *,'reclen needed exceeds integer 4-byte limit: ',b_reclen_potential
print *,' ',CUSTOM_REAL, NGLLSQUARE, num_abs_boundary_faces
print*,'bit size fortran: ',bit_size(b_reclen_potential)
@@ -862,7 +862,7 @@
filesize = filesize*NSTEP
! debug check size limit
- !if( NSTEP > 2147483647 / b_reclen_potential ) then
+ !if( NSTEP > 2147483646 / b_reclen_potential ) then
! print *,'file size needed exceeds integer 4-byte limit: ',b_reclen_potential,NSTEP
! print *,' ',CUSTOM_REAL, NGLLSQUARE, num_abs_boundary_faces,NSTEP
! print*,'file size fortran: ',filesize
@@ -909,7 +909,7 @@
! check integer size limit: size of b_reclen_field must fit onto an
! 4-byte integer
- if( num_abs_boundary_faces > 2147483647 / (CUSTOM_REAL * NDIM * NGLLSQUARE) ) then
+ if( num_abs_boundary_faces > 2147483646 / (CUSTOM_REAL * NDIM * NGLLSQUARE) ) then
print *,'reclen needed exceeds integer 4-byte limit: ',b_reclen_field_poro
print *,' ',CUSTOM_REAL, NDIM, NGLLSQUARE, num_abs_boundary_faces
print*,'bit size fortran: ',bit_size(b_reclen_field_poro)
Modified: seismo/3D/SPECFEM3D/trunk/src/specfem3D/save_external_bin_m_up.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/specfem3D/save_external_bin_m_up.f90 2012-10-13 18:42:30 UTC (rev 20836)
+++ seismo/3D/SPECFEM3D/trunk/src/specfem3D/save_external_bin_m_up.f90 2012-10-13 22:18:10 UTC (rev 20837)
@@ -73,7 +73,7 @@
use specfem_par_poroelastic,only: &
num_coupling_el_po_faces, &
nspec_inner_poroelastic,nspec_outer_poroelastic,num_phase_ispec_poroelastic,phase_ispec_inner_poroelastic
-
+
implicit none
include "constants.h"
@@ -229,7 +229,7 @@
write(IOUT) rmassz_acoustic
endif
endif
-
+
! free surface
write(IOUT) num_free_surface_faces
if( num_free_surface_faces > 0 ) then
Modified: seismo/3D/SPECFEM3D/trunk/src/specfem3D/setup_sources_receivers.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/specfem3D/setup_sources_receivers.f90 2012-10-13 18:42:30 UTC (rev 20836)
+++ seismo/3D/SPECFEM3D/trunk/src/specfem3D/setup_sources_receivers.f90 2012-10-13 22:18:10 UTC (rev 20837)
@@ -609,7 +609,7 @@
endif
! elastic source
if( ispec_is_elastic(ispec) ) then
- ! we use an inclined force defined by its magnitude and the projections
+ ! we use an inclined force defined by its magnitude and the projections
! of an arbitrary (non-unitary) direction vector on the E/N/Z_UP basis:
sourcearray(:,i,j,k) = factor_force_source(isource) * &
( nu_source(1,:,isource) * comp_dir_vect_source_E(isource) + &
@@ -878,7 +878,7 @@
etal = eta_source(isource)
gammal = gamma_source(isource)
endif
- call eval_shape3D_single(myrank,shape3D,xil,etal,gammal)
+ call eval_shape3D_single(myrank,shape3D,xil,etal,gammal,NGNOD)
! interpolates source locations
xmesh = 0.0
@@ -923,7 +923,7 @@
xil = xi_receiver(irec)
etal = eta_receiver(irec)
gammal = gamma_receiver(irec)
- call eval_shape3D_single(myrank,shape3D,xil,etal,gammal)
+ call eval_shape3D_single(myrank,shape3D,xil,etal,gammal,NGNOD)
! interpolates receiver locations
xmesh = 0.0
More information about the CIG-COMMITS
mailing list