[cig-commits] r20609 - in seismo/3D/SPECFEM3D/trunk/src: generate_databases meshfem3D shared specfem3D
dkomati1 at geodynamics.org
dkomati1 at geodynamics.org
Mon Aug 20 10:40:52 PDT 2012
Author: dkomati1
Date: 2012-08-20 10:40:52 -0700 (Mon, 20 Aug 2012)
New Revision: 20609
Modified:
seismo/3D/SPECFEM3D/trunk/src/generate_databases/get_coupling_surfaces.f90
seismo/3D/SPECFEM3D/trunk/src/meshfem3D/store_boundaries.f90
seismo/3D/SPECFEM3D/trunk/src/shared/constants.h.in
seismo/3D/SPECFEM3D/trunk/src/shared/detect_surface.f90
seismo/3D/SPECFEM3D/trunk/src/shared/get_element_face.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/recompute_jacobian.f90
seismo/3D/SPECFEM3D/trunk/src/specfem3D/setup_sources_receivers.f90
Log:
more modifications towards HEX27 support
Modified: seismo/3D/SPECFEM3D/trunk/src/generate_databases/get_coupling_surfaces.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/generate_databases/get_coupling_surfaces.f90 2012-08-20 17:12:38 UTC (rev 20608)
+++ seismo/3D/SPECFEM3D/trunk/src/generate_databases/get_coupling_surfaces.f90 2012-08-20 17:40:52 UTC (rev 20609)
@@ -173,7 +173,7 @@
! local parameters
! (assumes NGLLX=NGLLY=NGLLZ)
- 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(:,:,:),allocatable :: tmp_normal
@@ -182,7 +182,7 @@
integer,dimension(:,:,:),allocatable :: tmp_ijk
integer,dimension(:),allocatable :: tmp_ispec
- integer,dimension(NGNOD2D) :: iglob_corners_ref
+ integer,dimension(NGNOD2D_FOUR_CORNERS) :: iglob_corners_ref
integer :: ispec,i,j,igll,ier
integer :: inum,iface_ref
@@ -347,7 +347,7 @@
! local parameters
! (assumes NGLLX=NGLLY=NGLLZ)
- 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(:,:,:),allocatable :: tmp_normal
@@ -356,7 +356,7 @@
integer,dimension(:,:,:),allocatable :: tmp_ijk
integer,dimension(:),allocatable :: tmp_ispec
- integer,dimension(NGNOD2D) :: iglob_corners_ref
+ integer,dimension(NGNOD2D_FOUR_CORNERS) :: iglob_corners_ref
integer :: ispec,i,j,igll,ier
integer :: inum,iface_ref
@@ -512,7 +512,7 @@
! local parameters
! (assumes NGLLX=NGLLY=NGLLZ)
- 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(:,:,:),allocatable :: tmp_normal
@@ -521,7 +521,7 @@
integer,dimension(:,:,:),allocatable :: tmp_ijk,tmp_ijk_el
integer,dimension(:),allocatable :: tmp_ispec,tmp_ispec_el
- integer,dimension(NGNOD2D) :: iglob_corners_ref,iglob_corners_ref_el
+ integer,dimension(NGNOD2D_FOUR_CORNERS) :: iglob_corners_ref,iglob_corners_ref_el
integer :: ispec,i,j,k,igll,ier
integer :: ispec_el,ispec_ref_el
integer :: inum,iface_ref,iface_ref_el,iface_el,icorner
@@ -596,7 +596,7 @@
if(ispec_is_elastic(ispec_el))then
do iface_el=6,1,-1
! takes indices of corners of reference face
- do icorner = 1,NGNOD2D
+ do icorner = 1,NGNOD2D_FOUR_CORNERS
i = iface_all_corner_ijk(1,icorner,iface_el)
j = iface_all_corner_ijk(2,icorner,iface_el)
k = iface_all_corner_ijk(3,icorner,iface_el)
Modified: seismo/3D/SPECFEM3D/trunk/src/meshfem3D/store_boundaries.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/meshfem3D/store_boundaries.f90 2012-08-20 17:12:38 UTC (rev 20608)
+++ seismo/3D/SPECFEM3D/trunk/src/meshfem3D/store_boundaries.f90 2012-08-20 17:40:52 UTC (rev 20609)
@@ -50,7 +50,7 @@
! counters to keep track of number of elements on each of the boundaries
integer ispecb1,ispecb2,ispecb3,ispecb4,ispecb5,ispecb6
- ! check that the parameter is correct
+ ! check that the parameter file is correct
if(NGNOD2D_FOUR_CORNERS /= 4) call exit_MPI(myrank,'surface elements should have 4 control nodes in our internal mesher')
! initializes
@@ -183,7 +183,7 @@
double precision xelm(NGNOD2D_FOUR_CORNERS),yelm(NGNOD2D_FOUR_CORNERS),zelm(NGNOD2D_FOUR_CORNERS)
- ! check that the parameter is correct
+ ! check that the parameter file is correct
if(NGNOD2D_FOUR_CORNERS /= 4) call exit_MPI(myrank,'surface elements should have 4 control nodes in our internal mesher')
ispecb1 = 0
Modified: seismo/3D/SPECFEM3D/trunk/src/shared/constants.h.in
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/shared/constants.h.in 2012-08-20 17:12:38 UTC (rev 20608)
+++ seismo/3D/SPECFEM3D/trunk/src/shared/constants.h.in 2012-08-20 17:40:52 UTC (rev 20609)
@@ -304,7 +304,7 @@
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
+ double precision, parameter :: ZERO = 0.d0,ONE = 1.d0,TWO = 2.d0,HALF = 0.5d0
real(kind=CUSTOM_REAL), parameter :: &
ONE_THIRD = 1._CUSTOM_REAL/3._CUSTOM_REAL, &
Modified: seismo/3D/SPECFEM3D/trunk/src/shared/detect_surface.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/shared/detect_surface.f90 2012-08-20 17:12:38 UTC (rev 20608)
+++ seismo/3D/SPECFEM3D/trunk/src/shared/detect_surface.f90 2012-08-20 17:40:52 UTC (rev 20609)
@@ -260,8 +260,8 @@
real(kind=CUSTOM_REAL),dimension(6) :: midpoint_faces_x,midpoint_faces_y, &
midpoint_faces_z
real(kind=CUSTOM_REAL),dimension(6) :: midpoint_dist_x,midpoint_dist_y,midpoint_dist_z
- real(kind=CUSTOM_REAL),dimension(NGNOD2D) :: xcoord_face,ycoord_face,zcoord_face
- real(kind=CUSTOM_REAL) :: mindist,normal(NDIM)
+ real(kind=CUSTOM_REAL),dimension(NGNOD2D_FOUR_CORNERS) :: xcoord_face,ycoord_face,zcoord_face
+ real(kind=CUSTOM_REAL) :: min_dist,normal(NDIM)
integer, dimension(:), allocatable :: valence_external_mesh
integer :: ispec,i,j,k,iglob,ier,count
@@ -293,10 +293,10 @@
if( ier /= 0 ) stop 'error allocate valence array'
! an estimation of the minimum distance between global points (for an element width)
- mindist = minval( (xstore(ibool(1,3,3,:)) - xstore(ibool(NGLLX,3,3,:)))**2 &
+ min_dist = minval( (xstore(ibool(1,3,3,:)) - xstore(ibool(NGLLX,3,3,:)))**2 &
+ (ystore(ibool(1,3,3,:)) - ystore(ibool(NGLLX,3,3,:)))**2 &
+ (zstore(ibool(1,3,3,:)) - zstore(ibool(NGLLX,3,3,:)))**2 )
- mindist = sqrt(mindist)
+ min_dist = sqrt(min_dist)
! initialize surface indices
ispec_is_surface_external_mesh(:) = .false.
@@ -312,19 +312,19 @@
iglob = ibool(i,j,k,ispec)
! x cross-section
- if( abs( xstore(iglob) - x_section ) < 0.2*mindist ) then
+ if( abs( xstore(iglob) - x_section ) < 0.2*min_dist ) then
! sets valence to 1 for points on cross-sections
valence_external_mesh(iglob) = myrank+1
endif
! y cross-section
- if( abs( ystore(iglob) - y_section ) < 0.2*mindist ) then
+ if( abs( ystore(iglob) - y_section ) < 0.2*min_dist ) then
! sets valence to 1 for points on cross-sections
valence_external_mesh(iglob) = myrank+1
endif
! z cross-section
- if( abs( zstore(iglob) - z_section ) < 0.2*mindist ) then
+ if( abs( zstore(iglob) - z_section ) < 0.2*min_dist ) then
! sets valence to 1 for points on cross-sections
valence_external_mesh(iglob) = myrank+1
endif
@@ -391,17 +391,19 @@
if( ispec_has_points(ispec) ) then
! an estimation of the element width
- mindist = sqrt((xstore(ibool(1,3,3,ispec)) - xstore(ibool(NGLLX,3,3,ispec)))**2 &
- + (ystore(ibool(1,3,3,ispec)) - ystore(ibool(NGLLX,3,3,ispec)))**2 &
- + (zstore(ibool(1,3,3,ispec)) - zstore(ibool(NGLLX,3,3,ispec)))**2 )
+ min_dist = sqrt((xstore(ibool(1,3,3,ispec)) - xstore(ibool(NGLLX,3,3,ispec)))**2 &
+ + (ystore(ibool(1,3,3,ispec)) - ystore(ibool(NGLLX,3,3,ispec)))**2 &
+ + (zstore(ibool(1,3,3,ispec)) - zstore(ibool(NGLLX,3,3,ispec)))**2 )
! determines element face by minimum distance of midpoints
midpoint_faces_x(:) = 0.0
midpoint_faces_y(:) = 0.0
midpoint_faces_z(:) = 0.0
+
do iface=1,6
+
! face corners
- 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)
@@ -416,8 +418,8 @@
midpoint_faces_x(iface) = midpoint_faces_x(iface) + xcoord_face(icorner)
midpoint_faces_y(iface) = midpoint_faces_y(iface) + ycoord_face(icorner)
midpoint_faces_z(iface) = midpoint_faces_z(iface) + zcoord_face(icorner)
+ enddo
- enddo
midpoint_faces_x(iface) = midpoint_faces_x(iface) / 4.0
midpoint_faces_y(iface) = midpoint_faces_y(iface) / 4.0
midpoint_faces_z(iface) = midpoint_faces_z(iface) / 4.0
@@ -438,7 +440,7 @@
i = iface_midpoint_ijk(1,iface)
j = iface_midpoint_ijk(2,iface)
k = iface_midpoint_ijk(3,iface)
- if( midpoint_dist_x(iface) < 0.5*mindist .and. &
+ if( midpoint_dist_x(iface) < 0.5*min_dist .and. &
valence_external_mesh(ibool(i,j,k,ispec)) /= -1 ) then
! checks face normal points in similar direction as cross-section normal
if( abs(normal(1)) > 0.6 ) then
@@ -454,7 +456,7 @@
i = iface_midpoint_ijk(1,iface)
j = iface_midpoint_ijk(2,iface)
k = iface_midpoint_ijk(3,iface)
- if( midpoint_dist_y(iface) < 0.5*mindist .and. &
+ if( midpoint_dist_y(iface) < 0.5*min_dist .and. &
valence_external_mesh(ibool(i,j,k,ispec)) /= -1) then
! checks face normal points in similar direction as cross-section normal
if( abs(normal(2)) > 0.6 ) then
@@ -470,7 +472,7 @@
i = iface_midpoint_ijk(1,iface)
j = iface_midpoint_ijk(2,iface)
k = iface_midpoint_ijk(3,iface)
- if( midpoint_dist_z(iface) < 0.5*mindist .and. &
+ if( midpoint_dist_z(iface) < 0.5*min_dist .and. &
valence_external_mesh(ibool(i,j,k,ispec)) /= -1) then
! checks face normal points in similar direction as cross-section normal
if( abs(normal(3)) > 0.6 ) then
@@ -695,7 +697,7 @@
real(kind=CUSTOM_REAL), dimension(nglob) :: xstore,ystore,zstore
!local parameters
- real(kind=CUSTOM_REAL) :: mindist,distance
+ real(kind=CUSTOM_REAL) :: min_dist,distance
integer, dimension(:), allocatable :: valence_external_mesh
integer :: ispec,i,j,k,iglob,ier,count
real(kind=CUSTOM_REAL),parameter :: TOLERANCE_DISTANCE = 0.9
@@ -711,11 +713,11 @@
num_iglob_image_surface = 0
! an estimation of the minimum distance between global points
- mindist = minval( (xstore(ibool(1,1,1,:)) - xstore(ibool(2,1,1,:)))**2 &
- + (ystore(ibool(1,1,1,:)) - ystore(ibool(2,1,1,:)))**2 &
- + (zstore(ibool(1,1,1,:)) - zstore(ibool(2,1,1,:)))**2 )
- mindist = sqrt(mindist)
- distance = TOLERANCE_DISTANCE*mindist
+ min_dist = minval( (xstore(ibool(1,1,1,:)) - xstore(ibool(2,1,1,:)))**2 &
+ + (ystore(ibool(1,1,1,:)) - ystore(ibool(2,1,1,:)))**2 &
+ + (zstore(ibool(1,1,1,:)) - zstore(ibool(2,1,1,:)))**2 )
+ min_dist = sqrt(min_dist)
+ distance = TOLERANCE_DISTANCE*min_dist
! sets valence value to one corresponding to process rank for points on cross-sections
do ispec = 1, nspec
Modified: seismo/3D/SPECFEM3D/trunk/src/shared/get_element_face.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/shared/get_element_face.f90 2012-08-20 17:12:38 UTC (rev 20608)
+++ seismo/3D/SPECFEM3D/trunk/src/shared/get_element_face.f90 2012-08-20 17:40:52 UTC (rev 20609)
@@ -38,7 +38,7 @@
integer :: ispec,nspec,nglob,iface_id
! face corner locations
- real(kind=CUSTOM_REAL),dimension(NGNOD2D) :: xcoord,ycoord,zcoord
+ real(kind=CUSTOM_REAL),dimension(NGNOD2D_FOUR_CORNERS) :: xcoord,ycoord,zcoord
! index array
integer, dimension(NGLLX,NGLLY,NGLLZ,nspec) :: ibool
@@ -47,7 +47,7 @@
real(kind=CUSTOM_REAL) :: xstore_dummy(nglob),ystore_dummy(nglob),zstore_dummy(nglob)
! local parameters
- real(kind=CUSTOM_REAL),dimension(NGNOD2D) :: xcoord_face,ycoord_face,zcoord_face
+ real(kind=CUSTOM_REAL),dimension(NGNOD2D_FOUR_CORNERS) :: xcoord_face,ycoord_face,zcoord_face
real(kind=CUSTOM_REAL) :: midpoint_faces(NDIM,6),midpoint(NDIM),midpoint_distances(6)
! corners indices of reference cube faces
@@ -87,7 +87,7 @@
! gets face midpoint by its corners
midpoint(:) = 0.0
- do icorner=1,NGNOD2D
+ do icorner=1,NGNOD2D_FOUR_CORNERS
midpoint(1) = midpoint(1) + xcoord(icorner)
midpoint(2) = midpoint(2) + ycoord(icorner)
midpoint(3) = midpoint(3) + zcoord(icorner)
@@ -97,12 +97,12 @@
! determines element face by minimum distance of midpoints
midpoint_faces(:,:) = 0.0
do ifa=1,6
+
! face corners
- do icorner = 1,NGNOD2D
+ do icorner = 1,NGNOD2D_FOUR_CORNERS
i = iface_all_corner_ijk(1,icorner,ifa)
j = iface_all_corner_ijk(2,icorner,ifa)
k = iface_all_corner_ijk(3,icorner,ifa)
- !print*,'corner:',i,j,k,ispec
! coordinates
iglob = ibool(i,j,k,ispec)
@@ -115,6 +115,7 @@
midpoint_faces(2,ifa) = midpoint_faces(2,ifa) + ycoord_face(icorner)
midpoint_faces(3,ifa) = midpoint_faces(3,ifa) + zcoord_face(icorner)
enddo
+
midpoint_faces(:,ifa) = midpoint_faces(:,ifa) / 4.0
! distance
@@ -134,7 +135,7 @@
print*,'error element face midpoint distance:',midpoint_distances(iloc(1)), &
(xcoord(1)-xcoord(2))**2+(ycoord(1)-ycoord(2))**2+(zcoord(1)-zcoord(2))**2
! corner locations
- do icorner=1,NGNOD2D
+ do icorner=1,NGNOD2D_FOUR_CORNERS
i = iface_all_corner_ijk(1,icorner,iloc(1))
j = iface_all_corner_ijk(2,icorner,iloc(1))
k = iface_all_corner_ijk(3,icorner,iloc(1))
@@ -287,7 +288,7 @@
integer :: ispec,iface,nspec,nglob
! face corner locations
- real(kind=CUSTOM_REAL),dimension(NGNOD2D) :: xcoord,ycoord,zcoord
+ real(kind=CUSTOM_REAL),dimension(NGNOD2D_FOUR_CORNERS) :: xcoord,ycoord,zcoord
! index array
integer, dimension(NGLLX,NGLLY,NGLLZ,nspec) :: ibool
@@ -382,7 +383,7 @@
integer :: ispec,iface,nspec,nglob
! face corner locations
- real(kind=CUSTOM_REAL),dimension(NGNOD2D) :: xcoord,ycoord,zcoord
+ real(kind=CUSTOM_REAL),dimension(NGNOD2D_FOUR_CORNERS) :: xcoord,ycoord,zcoord
! index array
integer, dimension(NGLLX,NGLLY,NGLLZ,nspec) :: ibool
@@ -479,22 +480,22 @@
integer,intent(in) :: ispec,iface_ref,nspec,nglob
! face corner locations
- real(kind=CUSTOM_REAL),dimension(NGNOD2D),intent(out) :: xcoord,ycoord,zcoord
- integer,dimension(NGNOD2D),intent(out):: iglob_corners_ref
+ real(kind=CUSTOM_REAL),dimension(NGNOD2D_FOUR_CORNERS),intent(out) :: xcoord,ycoord,zcoord
+ integer,dimension(NGNOD2D_FOUR_CORNERS),intent(out):: iglob_corners_ref
! index array
integer, dimension(NGLLX,NGLLY,NGLLZ,nspec) :: ibool
! global point locations
real(kind=CUSTOM_REAL) :: xstore_dummy(nglob),ystore_dummy(nglob),zstore_dummy(nglob)
- ! assumes NGNOD2D == 4
+ ! assumes NGNOD2D_FOUR_CORNERS == 4
integer,dimension(3,4,6) :: iface_all_corner_ijk
! local parameters
integer :: icorner,i,j,k
! loops over corners
- do icorner = 1,NGNOD2D
+ do icorner = 1,NGNOD2D_FOUR_CORNERS
i = iface_all_corner_ijk(1,icorner,iface_ref)
j = iface_all_corner_ijk(2,icorner,iface_ref)
k = iface_all_corner_ijk(3,icorner,iface_ref)
Modified: seismo/3D/SPECFEM3D/trunk/src/shared/get_jacobian_boundaries.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/shared/get_jacobian_boundaries.f90 2012-08-20 17:12:38 UTC (rev 20608)
+++ seismo/3D/SPECFEM3D/trunk/src/shared/get_jacobian_boundaries.f90 2012-08-20 17:40:52 UTC (rev 20609)
@@ -258,152 +258,3 @@
end subroutine compute_jacobian_2D_face
-!
-!-------------------------
-!
-
-! This subroutine recompute the jacobian for one element
-! based upon the GLL points
-! Hejun Zhu, Oct 16, 2009
-
-! input: myrank,
-! xstore,ystore,zstore ----- input position
-! xigll,yigll,zigll ----- gll points position
-! ispec,nspec ----- element number
-! ACTUALLY_STORE_ARRAYS ------ save array or not
-
-! output: xixstore,xiystore,xizstore,
-! etaxstore,etaystore,etazstore,
-! gammaxstore,gammaystore,gammazstore ------ parameters used for calculating jacobian
-
-
- subroutine recalc_jacobian_gll2D(myrank,xstore,ystore,zstore, &
- xigll,yigll,wgllwgll,NGLLA,NGLLB, &
- ispec,nspec,jacobian2Dw_face,normal_face)
-
- implicit none
-
- include "constants.h"
-
- ! input parameter
- integer::myrank,ispec,nspec
- double precision, dimension(NGLLX,NGLLY,NGLLZ,nspec) :: xstore,ystore,zstore
-
- integer :: NGLLA,NGLLB
- double precision, dimension(NGLLA):: xigll
- double precision, dimension(NGLLB):: yigll
- double precision:: wgllwgll(NGLLA,NGLLB)
-
- real(kind=CUSTOM_REAL) jacobian2Dw_face(NGLLA,NGLLB)
- real(kind=CUSTOM_REAL) normal_face(NDIM,NGLLA,NGLLB)
-
- ! other parameters for this subroutine
- integer:: i,j,k,i1,j1,k1
- double precision:: xxi,xeta,yxi,yeta,zxi,zeta
- double precision:: xi,eta
- double precision,dimension(NGLLA):: hxir,hpxir
- double precision,dimension(NGLLB):: hetar,hpetar
- double precision:: hlagrange,hlagrange_xi,hlagrange_eta
- double precision:: jacobian
- double precision:: unx,uny,unz
-
- ! test parameters which can be deleted
- double precision:: xmesh,ymesh,zmesh
- double precision:: sumshape,sumdershapexi,sumdershapeeta
-
- ! first go over all gll points on face
- k=1
- do j=1,NGLLB
- do i=1,NGLLA
-
- xxi = 0.0
- xeta = 0.0
- yxi = 0.0
- yeta = 0.0
- zxi = 0.0
- zeta = 0.0
-
- xi = xigll(i)
- eta = yigll(j)
-
- ! calculate lagrange polynomial and its derivative
- call lagrange_any(xi,NGLLA,xigll,hxir,hpxir)
- call lagrange_any(eta,NGLLB,yigll,hetar,hpetar)
-
- ! test parameters
- sumshape = 0.0
- sumdershapexi = 0.0
- sumdershapeeta = 0.0
- xmesh = 0.0
- ymesh = 0.0
- zmesh = 0.0
-
- k1=1
- do j1 = 1,NGLLB
- do i1 = 1,NGLLA
- hlagrange = hxir(i1)*hetar(j1)
- hlagrange_xi = hpxir(i1)*hetar(j1)
- hlagrange_eta = hxir(i1)*hpetar(j1)
-
-
- xxi = xxi + xstore(i1,j1,k1,ispec)*hlagrange_xi
- xeta = xeta + xstore(i1,j1,k1,ispec)*hlagrange_eta
-
- yxi = yxi + ystore(i1,j1,k1,ispec)*hlagrange_xi
- yeta = yeta + ystore(i1,j1,k1,ispec)*hlagrange_eta
-
- zxi = zxi + zstore(i1,j1,k1,ispec)*hlagrange_xi
- zeta = zeta + zstore(i1,j1,k1,ispec)*hlagrange_eta
-
- ! test the lagrange polynomial and its derivate
- xmesh = xmesh + xstore(i1,j1,k1,ispec)*hlagrange
- ymesh = ymesh + ystore(i1,j1,k1,ispec)*hlagrange
- zmesh = zmesh + zstore(i1,j1,k1,ispec)*hlagrange
- sumshape = sumshape + hlagrange
- sumdershapexi = sumdershapexi + hlagrange_xi
- sumdershapeeta = sumdershapeeta + hlagrange_eta
-
- end do
- end do
-
- ! Check the lagrange polynomial and its derivative
- if (xmesh /=xstore(i,j,k,ispec).or.ymesh/=ystore(i,j,k,ispec).or.zmesh/=zstore(i,j,k,ispec)) then
- call exit_MPI(myrank,'new mesh positions are wrong in recalc_jacobian_gall3D.f90')
- end if
- if(abs(sumshape-one) > TINYVAL) then
- call exit_MPI(myrank,'error shape functions in recalc_jacobian_gll3D.f90')
- end if
- if(abs(sumdershapexi) > TINYVAL) then
- call exit_MPI(myrank,'error derivative xi shape functions in recalc_jacobian_gll3D.f90')
- end if
- if(abs(sumdershapeeta) > TINYVAL) then
- call exit_MPI(myrank,'error derivative eta shape functions in recalc_jacobian_gll3D.f90')
- end if
-
-! calculate the unnormalized normal to the boundary
- unx=yxi*zeta-yeta*zxi
- uny=zxi*xeta-zeta*xxi
- unz=xxi*yeta-xeta*yxi
- jacobian=dsqrt(unx**2+uny**2+unz**2)
- if(jacobian <= ZERO) call exit_MPI(myrank,'2D Jacobian undefined')
-
-! normalize normal vector and store weighted surface jacobian
-
-! distinguish if single or double precision for reals
- if(CUSTOM_REAL == SIZE_REAL) then
- 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)
- else
- jacobian2Dw_face(i,j) = jacobian * wgllwgll(i,j)
- normal_face(1,i,j)=unx/jacobian
- normal_face(2,i,j)=uny/jacobian
- normal_face(3,i,j)=unz/jacobian
- endif
-
- enddo
- enddo
-
- end subroutine recalc_jacobian_gll2D
-
Modified: seismo/3D/SPECFEM3D/trunk/src/shared/get_shape2D.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/shared/get_shape2D.f90 2012-08-20 17:12:38 UTC (rev 20608)
+++ seismo/3D/SPECFEM3D/trunk/src/shared/get_shape2D.f90 2012-08-20 17:40:52 UTC (rev 20609)
@@ -46,14 +46,18 @@
! location of the nodes of the 2D quadrilateral elements
double precision xi,eta
double precision xi_map,eta_map
+ double precision l1xi,l2xi,l3xi,l1eta,l2eta,l3eta
+ double precision l1pxi,l2pxi,l3pxi,l1peta,l2peta,l3peta
! for checking the 2D shape functions
double precision sumshape,sumdershapexi,sumdershapeeta
! 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(NGNOD /= 8 .and. NGNOD /= 27) call exit_MPI(myrank,'elements should have 8 or 27 control nodes')
+ if(NGNOD2D /= 4 .and. NGNOD2D /= 9) call exit_MPI(myrank,'surface elements should have 4 or 9 control nodes')
+ if(NGNOD2D == 4) then
+
! generate the 2D shape functions and their derivatives (4 nodes)
do i=1,NGLLA
@@ -88,6 +92,79 @@
enddo
enddo
+ else
+
+! generate the 2D shape functions and their derivatives (9 nodes)
+ do i=1,NGLLA
+
+ xi=xigll(i)
+
+ l1xi=HALF*xi*(xi-ONE)
+ l2xi=ONE-xi**2
+ l3xi=HALF*xi*(xi+ONE)
+
+ l1pxi=xi-HALF
+ l2pxi=-TWO*xi
+ l3pxi=xi+HALF
+
+ do j=1,NGLLB
+
+ eta=yigll(j)
+
+ l1eta=HALF*eta*(eta-ONE)
+ l2eta=ONE-eta**2
+ l3eta=HALF*eta*(eta+ONE)
+
+ l1peta=eta-HALF
+ l2peta=-TWO*eta
+ l3peta=eta+HALF
+
+! corner nodes
+
+ shape2D(1,i,j)=l1xi*l1eta
+ shape2D(2,i,j)=l3xi*l1eta
+ shape2D(3,i,j)=l3xi*l3eta
+ shape2D(4,i,j)=l1xi*l3eta
+
+ dershape2D(1,1,i,j)=l1pxi*l1eta
+ dershape2D(1,2,i,j)=l3pxi*l1eta
+ dershape2D(1,3,i,j)=l3pxi*l3eta
+ dershape2D(1,4,i,j)=l1pxi*l3eta
+
+ dershape2D(2,1,i,j)=l1xi*l1peta
+ dershape2D(2,2,i,j)=l3xi*l1peta
+ dershape2D(2,3,i,j)=l3xi*l3peta
+ dershape2D(2,4,i,j)=l1xi*l3peta
+
+! midside nodes
+
+ shape2D(5,i,j)=l2xi*l1eta
+ shape2D(6,i,j)=l3xi*l2eta
+ shape2D(7,i,j)=l2xi*l3eta
+ shape2D(8,i,j)=l1xi*l2eta
+
+ dershape2D(1,5,i,j)=l2pxi*l1eta
+ dershape2D(1,6,i,j)=l3pxi*l2eta
+ dershape2D(1,7,i,j)=l2pxi*l3eta
+ dershape2D(1,8,i,j)=l1pxi*l2eta
+
+ dershape2D(2,5,i,j)=l2xi*l1peta
+ dershape2D(2,6,i,j)=l3xi*l2peta
+ dershape2D(2,7,i,j)=l2xi*l3peta
+ dershape2D(2,8,i,j)=l1xi*l2peta
+
+! center node
+
+ shape2D(9,i,j)=l2xi*l2eta
+
+ dershape2D(1,9,i,j)=l2pxi*l2eta
+ dershape2D(2,9,i,j)=l2xi*l2peta
+
+ enddo
+ enddo
+
+ endif
+
! check the 2D shape functions
do i=1,NGLLA
do j=1,NGLLB
Modified: seismo/3D/SPECFEM3D/trunk/src/shared/get_shape3D.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/shared/get_shape3D.f90 2012-08-20 17:12:38 UTC (rev 20608)
+++ seismo/3D/SPECFEM3D/trunk/src/shared/get_shape3D.f90 2012-08-20 17:40:52 UTC (rev 20609)
@@ -157,7 +157,7 @@
! 3D shape functions for given, single xi/eta/gamma location
- subroutine get_shape3D_single(myrank,shape3D,xi,eta,gamma)
+ subroutine eval_shape3D_single(myrank,shape3D,xi,eta,gamma)
implicit none
@@ -210,13 +210,13 @@
! sum of derivative of shape functions should be zero
if(abs(sumshape-one) > TINYVAL) call exit_MPI(myrank,'error single shape functions')
- end subroutine get_shape3D_single
+ end subroutine eval_shape3D_single
!
!-------------------------------------------------------------------------------------------------
!
- subroutine get_shape3D_element_corners(xelm,yelm,zelm,ispec,&
+ subroutine eval_shape3D_element_corners(xelm,yelm,zelm,ispec,&
ibool,xstore,ystore,zstore,NSPEC_AB,NGLOB_AB)
implicit none
@@ -265,5 +265,5 @@
yelm(8)=ystore(ibool(1,NGLLY,NGLLZ,ispec))
zelm(8)=zstore(ibool(1,NGLLY,NGLLZ,ispec))
- end subroutine get_shape3D_element_corners
+ end subroutine eval_shape3D_element_corners
Modified: seismo/3D/SPECFEM3D/trunk/src/shared/recompute_jacobian.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/shared/recompute_jacobian.f90 2012-08-20 17:12:38 UTC (rev 20608)
+++ seismo/3D/SPECFEM3D/trunk/src/shared/recompute_jacobian.f90 2012-08-20 17:40:52 UTC (rev 20609)
@@ -48,6 +48,13 @@
double precision xgamma,ygamma,zgamma
double precision ra1,ra2,rb1,rb2,rc1,rc2
+ double precision l1xi,l2xi,l3xi
+ double precision l1eta,l2eta,l3eta
+ double precision l1gamma,l2gamma,l3gamma
+ double precision l1pxi,l2pxi,l3pxi
+ double precision l1peta,l2peta,l3peta
+ double precision l1pgamma,l2pgamma,l3pgamma
+
integer ia
! for 8-node element
@@ -56,8 +63,10 @@
! recompute jacobian for any (xi,eta,gamma) point, not necessarily a GLL point
! check that the parameter file is correct
- if(NGNOD /= 8) stop 'elements should have 8 control nodes'
+ if(NGNOD /= 8 .and. NGNOD /= 27) stop 'elements should have 8 or 27 control nodes'
+ if(NGNOD == 8) then
+
! ***
! *** create the 3D shape functions and the Jacobian for an 8-node element
! ***
@@ -109,6 +118,168 @@
dershape3D(3,7) = ONE_EIGHTH*ra1*rb1
dershape3D(3,8) = ONE_EIGHTH*ra2*rb1
+ else
+
+! ***
+! *** create the 3D shape functions and the Jacobian for a 27-node element
+! ***
+
+ 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)=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
+
+ dershape3D(1,1)=l1pxi*l1eta*l1gamma
+ dershape3D(1,2)=l3pxi*l1eta*l1gamma
+ dershape3D(1,3)=l3pxi*l3eta*l1gamma
+ dershape3D(1,4)=l1pxi*l3eta*l1gamma
+ dershape3D(1,5)=l1pxi*l1eta*l3gamma
+ dershape3D(1,6)=l3pxi*l1eta*l3gamma
+ dershape3D(1,7)=l3pxi*l3eta*l3gamma
+ dershape3D(1,8)=l1pxi*l3eta*l3gamma
+
+ dershape3D(2,1)=l1xi*l1peta*l1gamma
+ dershape3D(2,2)=l3xi*l1peta*l1gamma
+ dershape3D(2,3)=l3xi*l3peta*l1gamma
+ dershape3D(2,4)=l1xi*l3peta*l1gamma
+ dershape3D(2,5)=l1xi*l1peta*l3gamma
+ dershape3D(2,6)=l3xi*l1peta*l3gamma
+ dershape3D(2,7)=l3xi*l3peta*l3gamma
+ dershape3D(2,8)=l1xi*l3peta*l3gamma
+
+ dershape3D(3,1)=l1xi*l1eta*l1pgamma
+ dershape3D(3,2)=l3xi*l1eta*l1pgamma
+ dershape3D(3,3)=l3xi*l3eta*l1pgamma
+ dershape3D(3,4)=l1xi*l3eta*l1pgamma
+ dershape3D(3,5)=l1xi*l1eta*l3pgamma
+ dershape3D(3,6)=l3xi*l1eta*l3pgamma
+ dershape3D(3,7)=l3xi*l3eta*l3pgamma
+ dershape3D(3,8)=l1xi*l3eta*l3pgamma
+
+! 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
+
+ dershape3D(1,9)=l2pxi*l1eta*l1gamma
+ dershape3D(1,10)=l3pxi*l2eta*l1gamma
+ dershape3D(1,11)=l2pxi*l3eta*l1gamma
+ dershape3D(1,12)=l1pxi*l2eta*l1gamma
+ dershape3D(1,13)=l1pxi*l1eta*l2gamma
+ dershape3D(1,14)=l3pxi*l1eta*l2gamma
+ dershape3D(1,15)=l3pxi*l3eta*l2gamma
+ dershape3D(1,16)=l1pxi*l3eta*l2gamma
+ dershape3D(1,17)=l2pxi*l1eta*l3gamma
+ dershape3D(1,18)=l3pxi*l2eta*l3gamma
+ dershape3D(1,19)=l2pxi*l3eta*l3gamma
+ dershape3D(1,20)=l1pxi*l2eta*l3gamma
+
+ dershape3D(2,9)=l2xi*l1peta*l1gamma
+ dershape3D(2,10)=l3xi*l2peta*l1gamma
+ dershape3D(2,11)=l2xi*l3peta*l1gamma
+ dershape3D(2,12)=l1xi*l2peta*l1gamma
+ dershape3D(2,13)=l1xi*l1peta*l2gamma
+ dershape3D(2,14)=l3xi*l1peta*l2gamma
+ dershape3D(2,15)=l3xi*l3peta*l2gamma
+ dershape3D(2,16)=l1xi*l3peta*l2gamma
+ dershape3D(2,17)=l2xi*l1peta*l3gamma
+ dershape3D(2,18)=l3xi*l2peta*l3gamma
+ dershape3D(2,19)=l2xi*l3peta*l3gamma
+ dershape3D(2,20)=l1xi*l2peta*l3gamma
+
+ dershape3D(3,9)=l2xi*l1eta*l1pgamma
+ dershape3D(3,10)=l3xi*l2eta*l1pgamma
+ dershape3D(3,11)=l2xi*l3eta*l1pgamma
+ dershape3D(3,12)=l1xi*l2eta*l1pgamma
+ dershape3D(3,13)=l1xi*l1eta*l2pgamma
+ dershape3D(3,14)=l3xi*l1eta*l2pgamma
+ dershape3D(3,15)=l3xi*l3eta*l2pgamma
+ dershape3D(3,16)=l1xi*l3eta*l2pgamma
+ dershape3D(3,17)=l2xi*l1eta*l3pgamma
+ dershape3D(3,18)=l3xi*l2eta*l3pgamma
+ dershape3D(3,19)=l2xi*l3eta*l3pgamma
+ dershape3D(3,20)=l1xi*l2eta*l3pgamma
+
+! 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
+
+ dershape3D(1,21)=l2pxi*l2eta*l1gamma
+ dershape3D(1,22)=l2pxi*l1eta*l2gamma
+ dershape3D(1,23)=l3pxi*l2eta*l2gamma
+ dershape3D(1,24)=l2pxi*l3eta*l2gamma
+ dershape3D(1,25)=l1pxi*l2eta*l2gamma
+ dershape3D(1,26)=l2pxi*l2eta*l3gamma
+
+ dershape3D(2,21)=l2xi*l2peta*l1gamma
+ dershape3D(2,22)=l2xi*l1peta*l2gamma
+ dershape3D(2,23)=l3xi*l2peta*l2gamma
+ dershape3D(2,24)=l2xi*l3peta*l2gamma
+ dershape3D(2,25)=l1xi*l2peta*l2gamma
+ dershape3D(2,26)=l2xi*l2peta*l3gamma
+
+ dershape3D(3,21)=l2xi*l2eta*l1pgamma
+ dershape3D(3,22)=l2xi*l1eta*l2pgamma
+ dershape3D(3,23)=l3xi*l2eta*l2pgamma
+ dershape3D(3,24)=l2xi*l3eta*l2pgamma
+ dershape3D(3,25)=l1xi*l2eta*l2pgamma
+ dershape3D(3,26)=l2xi*l2eta*l3pgamma
+
+! center node
+
+ shape3D(27)=l2xi*l2eta*l2gamma
+
+ dershape3D(1,27)=l2pxi*l2eta*l2gamma
+ dershape3D(2,27)=l2xi*l2peta*l2gamma
+ dershape3D(3,27)=l2xi*l2eta*l2pgamma
+
+ endif
+
! compute coordinates and jacobian matrix
x=ZERO
y=ZERO
Modified: seismo/3D/SPECFEM3D/trunk/src/specfem3D/setup_sources_receivers.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/specfem3D/setup_sources_receivers.f90 2012-08-20 17:12:38 UTC (rev 20608)
+++ seismo/3D/SPECFEM3D/trunk/src/specfem3D/setup_sources_receivers.f90 2012-08-20 17:40:52 UTC (rev 20609)
@@ -829,7 +829,7 @@
! gets element ancor nodes
if( myrank == islice_selected_source(isource) ) then
! find the coordinates of the eight corner nodes of the element
- call get_shape3D_element_corners(xelm,yelm,zelm,ispec,&
+ call eval_shape3D_element_corners(xelm,yelm,zelm,ispec,&
ibool,xstore,ystore,zstore,NSPEC_AB,NGLOB_AB)
endif
@@ -859,7 +859,7 @@
etal = eta_source(isource)
gammal = gamma_source(isource)
endif
- call get_shape3D_single(myrank,shape3D,xil,etal,gammal)
+ call eval_shape3D_single(myrank,shape3D,xil,etal,gammal)
! interpolates source locations
xmesh = 0.0
@@ -883,7 +883,7 @@
! find the coordinates of the eight corner nodes of the element
if( myrank == islice_selected_rec(irec) ) then
- call get_shape3D_element_corners(xelm,yelm,zelm,ispec,&
+ call eval_shape3D_element_corners(xelm,yelm,zelm,ispec,&
ibool,xstore,ystore,zstore,NSPEC_AB,NGLOB_AB)
endif
! master collects corner locations
@@ -904,7 +904,7 @@
xil = xi_receiver(irec)
etal = eta_receiver(irec)
gammal = gamma_receiver(irec)
- call get_shape3D_single(myrank,shape3D,xil,etal,gammal)
+ call eval_shape3D_single(myrank,shape3D,xil,etal,gammal)
! interpolates receiver locations
xmesh = 0.0
More information about the CIG-COMMITS
mailing list