[cig-commits] r20587 - in seismo/3D/SPECFEM3D/trunk/src: check_mesh_quality_CUBIT_Abaqus decompose_mesh_SCOTCH meshfem3D shared
dkomati1 at geodynamics.org
dkomati1 at geodynamics.org
Thu Aug 16 17:28:34 PDT 2012
Author: dkomati1
Date: 2012-08-16 17:28:33 -0700 (Thu, 16 Aug 2012)
New Revision: 20587
Modified:
seismo/3D/SPECFEM3D/trunk/src/check_mesh_quality_CUBIT_Abaqus/check_mesh_quality_CUBIT_Abaqus.f90
seismo/3D/SPECFEM3D/trunk/src/decompose_mesh_SCOTCH/part_decompose_mesh_SCOTCH.f90
seismo/3D/SPECFEM3D/trunk/src/meshfem3D/store_coords.f90
seismo/3D/SPECFEM3D/trunk/src/shared/create_movie_shakemap_AVS_DX_GMT.f90
seismo/3D/SPECFEM3D/trunk/src/shared/get_element_face.f90
seismo/3D/SPECFEM3D/trunk/src/shared/get_jacobian_boundaries.f90
Log:
removed obsolete comments and obsolete / unused code
Modified: seismo/3D/SPECFEM3D/trunk/src/check_mesh_quality_CUBIT_Abaqus/check_mesh_quality_CUBIT_Abaqus.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/check_mesh_quality_CUBIT_Abaqus/check_mesh_quality_CUBIT_Abaqus.f90 2012-08-17 00:10:01 UTC (rev 20586)
+++ seismo/3D/SPECFEM3D/trunk/src/check_mesh_quality_CUBIT_Abaqus/check_mesh_quality_CUBIT_Abaqus.f90 2012-08-17 00:28:33 UTC (rev 20587)
@@ -50,60 +50,6 @@
! number of points and of hex or quad elements
! number of points of a hex or quad element
-! character(len=100), parameter :: cubit_mesh_file = 'HOMOGENE_3D_lisse_300_in_meters.inp'
-! integer, parameter :: NPOIN = 98692, NSPEC = 90585, NGNOD = 8
-! logical, parameter :: IGNORE_OTHER_HEADERS = .false.
-! double precision, parameter :: delta_t = 1.d-3
-! double precision, parameter :: VP_MAX = 3000.d0
-
-! character(len=100), parameter :: cubit_mesh_file = 'regolite_3D_rego3d_70m_in_meters.inp'
-! integer, parameter :: NPOIN = 4050696, NSPEC = 3410265, NGNOD = 8
-! logical, parameter :: IGNORE_OTHER_HEADERS = .false.
-! double precision, parameter :: delta_t = 3.d-4
-! double precision, parameter :: VP_MAX = 900.d0 ! because the smallest element is in the regolith layer, not in the bedrock
-
-! character(len=100), parameter :: cubit_mesh_file = 'rego3d_70_disp.inp'
-! integer, parameter :: NPOIN = 5924713, NSPEC = 5797440, NGNOD = 8
-! logical, parameter :: IGNORE_OTHER_HEADERS = .true.
-! double precision, parameter :: delta_t = 3.d-4
-! double precision, parameter :: VP_MAX = 3000.d0
-
-! character(len=100), parameter :: cubit_mesh_file = 'rego3d_70_disp_regolith_only.inp'
-! integer, parameter :: NPOIN = 5924713, NSPEC = 252928, NGNOD = 8
-! logical, parameter :: IGNORE_OTHER_HEADERS = .false.
-! double precision, parameter :: delta_t = 3.d-4
-! double precision, parameter :: VP_MAX = 900.d0 ! because only regolith, no bedrock
-
-! character(len=100), parameter :: cubit_mesh_file = 'rego3d_70_disp_bedrock_only.inp'
-! integer, parameter :: NPOIN = 5924713, NSPEC = 5797440 - 252928, NGNOD = 8
-! logical, parameter :: IGNORE_OTHER_HEADERS = .false.
-! double precision, parameter :: delta_t = 3.d-4
-! double precision, parameter :: VP_MAX = 3000.d0
-
-! character(len=100), parameter :: cubit_mesh_file = 'HOMOGENE_2D_in_meters.inp'
-! integer, parameter :: NPOIN = 3882, NSPEC = 3744, NGNOD = 4
-! logical, parameter :: IGNORE_OTHER_HEADERS = .false.
-! double precision, parameter :: delta_t = 5.d-3
-! double precision, parameter :: VP_MAX = 3000.d0
-
-! character(len=100), parameter :: cubit_mesh_file = 'eros_complexe_2d_regolite_fractures_modifie_in_meters.inp'
-! integer, parameter :: NPOIN = 57807, NSPEC = 56983, NGNOD = 4
-! logical, parameter :: IGNORE_OTHER_HEADERS = .true.
-! double precision, parameter :: delta_t = 1.5d-4
-! double precision, parameter :: VP_MAX = 3000.d0
-
-! character(len=100), parameter :: cubit_mesh_file = 'REGOLITE_only_no_fractures_2D_in_meters.inp'
-! integer, parameter :: NPOIN = 32536, NSPEC = 31695, NGNOD = 4
-! logical, parameter :: IGNORE_OTHER_HEADERS = .true.
-! double precision, parameter :: delta_t = 1.5d-4
-! double precision, parameter :: VP_MAX = 900.d0 ! because the smallest element is in the regolith layer, not in the bedrock
-
-! character(len=100), parameter :: cubit_mesh_file = 'david_mesh_doubl_500_900_6layers.inp'
-! integer, parameter :: NPOIN = 4513255, NSPEC = 4379190, NGNOD = 8
-! logical, parameter :: IGNORE_OTHER_HEADERS = .false.
-! double precision, parameter :: delta_t = 1.5d-4
-! double precision, parameter :: VP_MAX = 3000.d0
-
! example: layered_halfspace
! Cubit -> File -> Export... Abacus (*.inp)
! ( block ids: 1 2 3 ) volumes only
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-08-17 00:10:01 UTC (rev 20586)
+++ seismo/3D/SPECFEM3D/trunk/src/decompose_mesh_SCOTCH/part_decompose_mesh_SCOTCH.f90 2012-08-17 00:28:33 UTC (rev 20587)
@@ -43,8 +43,6 @@
integer, parameter :: ELASTIC_LOAD = 4
integer, parameter :: POROELASTIC_LOAD = 8
-! include './constants_decompose_mesh_SCOTCH.h'
-
contains
!-----------------------------------------------
@@ -198,8 +196,6 @@
subroutine build_glob2loc_nodes(nspec, nnodes, nsize, nnodes_elmnts, nodes_elmnts, part, &
glob2loc_nodes_nparts, glob2loc_nodes_parts, glob2loc_nodes,nparts)
-! include './constants_decompose_mesh_SCOTCH.h'
-
integer, intent(in) :: nspec
integer, intent(in) :: nsize
integer, intent(in) :: nnodes
@@ -590,9 +586,6 @@
do i = 0, nnodes-1
do j = glob2loc_nodes_nparts(i), glob2loc_nodes_nparts(i+1)-1
if ( glob2loc_nodes_parts(j) == iproc ) then
- !write(IIN_database,*) glob2loc_nodes(j)+1, nodes_coords(1,i+1), &
- ! nodes_coords(2,i+1), nodes_coords(3,i+1)
-
write(IIN_database) glob2loc_nodes(j)+1, nodes_coords(1,i+1), &
nodes_coords(2,i+1), nodes_coords(3,i+1)
end if
@@ -615,8 +608,8 @@
character (len=30), dimension(6,count_undef_mat) :: undef_mat_prop
integer :: i
- !write(IIN_database,*) count_def_mat,count_undef_mat
write(IIN_database) count_def_mat,count_undef_mat
+
do i = 1, count_def_mat
! database material definition
!
@@ -626,24 +619,14 @@
!
! (note that this order of the properties is different than the input in nummaterial_velocity_file)
!
- !write(IIN_database,*) mat_prop(1,i), mat_prop(2,i), mat_prop(3,i), &
- ! mat_prop(4,i), mat_prop(5,i), mat_prop(6,i), &
- ! mat_prop(7,i), mat_prop(8,i), mat_prop(9,i), &
- ! mat_prop(10,i), mat_prop(11,i), mat_prop(12,i), &
- ! mat_prop(13,i), mat_prop(14,i), mat_prop(15,i), mat_prop(16,i)
-
write(IIN_database) mat_prop(1,i), mat_prop(2,i), mat_prop(3,i), &
mat_prop(4,i), mat_prop(5,i), mat_prop(6,i), &
mat_prop(7,i), mat_prop(8,i), mat_prop(9,i), &
mat_prop(10,i), mat_prop(11,i), mat_prop(12,i), &
mat_prop(13,i), mat_prop(14,i), mat_prop(15,i), mat_prop(16,i)
-
end do
- do i = 1, count_undef_mat
- !write(IIN_database,*) trim(undef_mat_prop(1,i)),' ',trim(undef_mat_prop(2,i)),' ', &
- ! trim(undef_mat_prop(3,i)),' ',trim(undef_mat_prop(4,i)),' ', &
- ! trim(undef_mat_prop(5,i)),' ',trim(undef_mat_prop(6,i))
+ do i = 1, count_undef_mat
write(IIN_database) undef_mat_prop(1,i),undef_mat_prop(2,i), &
undef_mat_prop(3,i),undef_mat_prop(4,i), &
undef_mat_prop(5,i),undef_mat_prop(6,i)
@@ -702,7 +685,6 @@
loc_nspec2D_xmin = loc_nspec2D_xmin + 1
end if
end do
- !write(IIN_database,*) 1, loc_nspec2D_xmin
write(IIN_database) 1, loc_nspec2D_xmin
loc_nspec2D_xmax = 0
@@ -711,7 +693,6 @@
loc_nspec2D_xmax = loc_nspec2D_xmax + 1
end if
end do
- !write(IIN_database,*) 2, loc_nspec2D_xmax
write(IIN_database) 2, loc_nspec2D_xmax
loc_nspec2D_ymin = 0
@@ -720,7 +701,6 @@
loc_nspec2D_ymin = loc_nspec2D_ymin + 1
end if
end do
- !write(IIN_database,*) 3, loc_nspec2D_ymin
write(IIN_database) 3, loc_nspec2D_ymin
loc_nspec2D_ymax = 0
@@ -729,7 +709,6 @@
loc_nspec2D_ymax = loc_nspec2D_ymax + 1
end if
end do
- !write(IIN_database,*) 4, loc_nspec2D_ymax
write(IIN_database) 4, loc_nspec2D_ymax
loc_nspec2D_bottom = 0
@@ -738,7 +717,6 @@
loc_nspec2D_bottom = loc_nspec2D_bottom + 1
end if
end do
- !write(IIN_database,*) 5, loc_nspec2D_bottom
write(IIN_database) 5, loc_nspec2D_bottom
loc_nspec2D_top = 0
@@ -747,7 +725,6 @@
loc_nspec2D_top = loc_nspec2D_top + 1
end if
end do
- !write(IIN_database,*) 6, loc_nspec2D_top
write(IIN_database) 6, loc_nspec2D_top
! outputs element index and element node indices
@@ -781,9 +758,6 @@
loc_node4 = glob2loc_nodes(j)+1
end if
end do
- !write(IIN_database,*) glob2loc_elmnts(ibelm_xmin(i)-1)+1, &
- ! loc_node1, loc_node2, loc_node3, loc_node4
-
write(IIN_database) glob2loc_elmnts(ibelm_xmin(i)-1)+1, &
loc_node1, loc_node2, loc_node3, loc_node4
end if
@@ -815,8 +789,6 @@
loc_node4 = glob2loc_nodes(j)+1
end if
end do
- !write(IIN_database,*) glob2loc_elmnts(ibelm_xmax(i)-1)+1, &
- ! loc_node1, loc_node2, loc_node3, loc_node4
write(IIN_database) glob2loc_elmnts(ibelm_xmax(i)-1)+1, &
loc_node1, loc_node2, loc_node3, loc_node4
end if
@@ -848,8 +820,6 @@
loc_node4 = glob2loc_nodes(j)+1
end if
end do
- !write(IIN_database,*) glob2loc_elmnts(ibelm_ymin(i)-1)+1, &
- ! loc_node1, loc_node2, loc_node3, loc_node4
write(IIN_database) glob2loc_elmnts(ibelm_ymin(i)-1)+1, &
loc_node1, loc_node2, loc_node3, loc_node4
end if
@@ -881,8 +851,6 @@
loc_node4 = glob2loc_nodes(j)+1
end if
end do
- !write(IIN_database,*) glob2loc_elmnts(ibelm_ymax(i)-1)+1, &
- ! loc_node1, loc_node2, loc_node3, loc_node4
write(IIN_database) glob2loc_elmnts(ibelm_ymax(i)-1)+1, &
loc_node1, loc_node2, loc_node3, loc_node4
end if
@@ -914,7 +882,6 @@
loc_node4 = glob2loc_nodes(j)+1
end if
end do
- !write(IIN_database,*) glob2loc_elmnts(ibelm_bottom(i)-1)+1, loc_node1, loc_node2, loc_node3, loc_node4
write(IIN_database) glob2loc_elmnts(ibelm_bottom(i)-1)+1, loc_node1, loc_node2, loc_node3, loc_node4
end if
end do
@@ -941,7 +908,6 @@
loc_node4 = glob2loc_nodes(j)+1
end if
end do
- !write(IIN_database,*) glob2loc_elmnts(ibelm_top(i)-1)+1, loc_node1, loc_node2, loc_node3, loc_node4
write(IIN_database) glob2loc_elmnts(ibelm_top(i)-1)+1, loc_node1, loc_node2, loc_node3, loc_node4
end if
@@ -1000,8 +966,6 @@
! format:
! # ispec_local # material_index_1 # material_index_2 # corner_id1 # corner_id2 # ... # corner_id8
- !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)
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
@@ -1022,8 +986,6 @@
glob2loc_nodes_nparts, glob2loc_nodes_parts, &
glob2loc_nodes, num_phase, nparts)
-! include './constants_decompose_mesh_SCOTCH.h'
-
integer, intent(in) :: IIN_database
integer, intent(in) :: iproc
integer, intent(in) :: ninterfaces, nparts
@@ -1076,10 +1038,8 @@
do j = i+1, nparts-1
if ( my_interfaces(num_interface) == 1 ) then
if ( i == iproc ) then
- !write(IIN_database,*) j, my_nb_interfaces(num_interface)
write(IIN_database) j, my_nb_interfaces(num_interface)
else
- !write(IIN_database,*) i, my_nb_interfaces(num_interface)
write(IIN_database) i, my_nb_interfaces(num_interface)
end if
@@ -1091,34 +1051,6 @@
local_elmnt = glob2loc_elmnts(tab_interfaces(k*7+1))+1
end if
-!!$ if ( tab_interfaces(k*7+2) == 1 ) then
-!!$ do l = glob2loc_nodes_nparts(tab_interfaces(k*7+3)), &
-!!$ glob2loc_nodes_nparts(tab_interfaces(k*7+3)+1)-1
-!!$ if ( glob2loc_nodes_parts(l) == iproc ) then
-!!$ local_nodes(1) = glob2loc_nodes(l)+1
-!!$ end if
-!!$ end do
-!!$
-!!$ write(IIN_database,*) local_elmnt, tab_interfaces(k*7+2), local_nodes(1), -1
-!!$ else
-!!$ if ( tab_interfaces(k*7+2) == 2 ) then
-!!$ do l = glob2loc_nodes_nparts(tab_interfaces(k*7+3)), &
-!!$ glob2loc_nodes_nparts(tab_interfaces(k*7+3)+1)-1
-!!$ if ( glob2loc_nodes_parts(l) == iproc ) then
-!!$ local_nodes(1) = glob2loc_nodes(l)+1
-!!$ end if
-!!$ end do
-!!$ do l = glob2loc_nodes_nparts(tab_interfaces(k*7+4)), &
-!!$ glob2loc_nodes_nparts(tab_interfaces(k*7+4)+1)-1
-!!$ if ( glob2loc_nodes_parts(l) == iproc ) then
-!!$ local_nodes(2) = glob2loc_nodes(l)+1
-!!$ end if
-!!$ end do
-!!$ write(IIN_database,*) local_elmnt, tab_interfaces(k*7+2), local_nodes(1), local_nodes(2)
-!!$ else
-!!$ write(IIN_database,*) "erreur_write_interface_", tab_interfaces(k*7+2)
-!!$ end if
-!!$ end if
select case (tab_interfaces(k*7+2))
case (1)
! single point element
@@ -1128,8 +1060,6 @@
local_nodes(1) = glob2loc_nodes(l)+1
end if
end do
- !write(IIN_database,*) local_elmnt, tab_interfaces(k*7+2), &
- ! local_nodes(1), -1, -1, -1
write(IIN_database) local_elmnt, tab_interfaces(k*7+2), &
local_nodes(1), -1, -1, -1
@@ -1147,8 +1077,6 @@
local_nodes(2) = glob2loc_nodes(l)+1
end if
end do
- !write(IIN_database,*) local_elmnt, tab_interfaces(k*7+2), &
- ! local_nodes(1), local_nodes(2), -1, -1
write(IIN_database) local_elmnt, tab_interfaces(k*7+2), &
local_nodes(1), local_nodes(2), -1, -1
@@ -1179,8 +1107,6 @@
local_nodes(4) = glob2loc_nodes(l)+1
end if
end do
- !write(IIN_database,*) local_elmnt, tab_interfaces(k*7+2), &
- ! local_nodes(1), local_nodes(2),local_nodes(3), local_nodes(4)
write(IIN_database) local_elmnt, tab_interfaces(k*7+2), &
local_nodes(1), local_nodes(2),local_nodes(3), local_nodes(4)
@@ -1242,7 +1168,6 @@
if( loc_nspec2D_moho == 0 ) return
! format: #surface_id, #number of elements
- !write(IIN_database,*) 7, loc_nspec2D_moho
write(IIN_database) 7, loc_nspec2D_moho
! outputs element index and element node indices
@@ -1274,8 +1199,6 @@
loc_node4 = glob2loc_nodes(j)+1
end if
end do
- !write(IIN_database,*) glob2loc_elmnts(ibelm_moho(i)-1)+1, &
- ! loc_node1, loc_node2, loc_node3, loc_node4
write(IIN_database) glob2loc_elmnts(ibelm_moho(i)-1)+1, &
loc_node1, loc_node2, loc_node3, loc_node4
end if
@@ -1436,7 +1359,6 @@
if( ier /= 0 ) stop 'error allocating array nnodes_elmnts'
allocate(nodes_elmnts(0:nsize*nnodes-1),stat=ier)
if( ier /= 0 ) stop 'error allocating array nodes_elmnts'
- !call mesh2dual_ncommonnodes(nspec, nnodes, elmnts, xadj, adjncy, nnodes_elmnts, nodes_elmnts,4)
call mesh2dual_ncommonnodes(nspec, nnodes, nsize, sup_neighbour, &
elmnts, xadj, adjncy, nnodes_elmnts, &
nodes_elmnts, max_neighbour, 4)
@@ -1664,6 +1586,5 @@
end subroutine moho_surface_repartitioning
-
end module part_decompose_mesh_SCOTCH
Modified: seismo/3D/SPECFEM3D/trunk/src/meshfem3D/store_coords.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/meshfem3D/store_coords.f90 2012-08-17 00:10:01 UTC (rev 20586)
+++ seismo/3D/SPECFEM3D/trunk/src/meshfem3D/store_coords.f90 2012-08-17 00:28:33 UTC (rev 20587)
@@ -32,113 +32,10 @@
integer ispec,nspec
-! double precision shape3D(NGNOD,NGLLX,NGLLY,NGLLZ)
-! double precision dershape3D(NDIM,NGNOD,NGLLX,NGLLY,NGLLZ)
-
double precision, dimension(NGNOD) :: xelm,yelm,zelm
-! real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,nspec) :: &
-! xixstore,xiystore,xizstore,etaxstore,etaystore,etazstore, &
-! gammaxstore,gammaystore,gammazstore,jacobianstore
-
double precision, dimension(NGLLX,NGLLY,NGLLZ,nspec) :: xstore,ystore,zstore
-! integer i,j,k,ia
-! double precision xxi,xeta,xgamma,yxi,yeta,ygamma,zxi,zeta,zgamma
-! double precision xmesh,ymesh,zmesh
-! double precision xix,xiy,xiz,etax,etay,etaz,gammax,gammay,gammaz
-! double precision jacobian
-
- ! do k=1,NGLLZ
-! do j=1,NGLLY
-! do i=1,NGLLX
-
-! ! xxi = ZERO
-! ! xeta = ZERO
-! ! xgamma = ZERO
-! ! yxi = ZERO
-! ! yeta = ZERO
-! ! ygamma = ZERO
-! ! zxi = ZERO
-! ! zeta = ZERO
-! ! zgamma = ZERO
-! xmesh = ZERO
-! ymesh = ZERO
-! zmesh = ZERO
-
-! do ia=1,NGNOD
-! ! xxi = xxi + dershape3D(1,ia,i,j,k)*xelm(ia)
-! ! xeta = xeta + dershape3D(2,ia,i,j,k)*xelm(ia)
-! ! xgamma = xgamma + dershape3D(3,ia,i,j,k)*xelm(ia)
-! ! yxi = yxi + dershape3D(1,ia,i,j,k)*yelm(ia)
-! ! yeta = yeta + dershape3D(2,ia,i,j,k)*yelm(ia)
-! ! ygamma = ygamma + dershape3D(3,ia,i,j,k)*yelm(ia)
-! ! zxi = zxi + dershape3D(1,ia,i,j,k)*zelm(ia)
-! ! zeta = zeta + dershape3D(2,ia,i,j,k)*zelm(ia)
-! ! zgamma = zgamma + dershape3D(3,ia,i,j,k)*zelm(ia)
-! xmesh = xmesh + shape3D(ia,i,j,k)*xelm(ia)
-! ymesh = ymesh + shape3D(ia,i,j,k)*yelm(ia)
-! zmesh = zmesh + shape3D(ia,i,j,k)*zelm(ia)
-! enddo
-
-! ! jacobian = xxi*(yeta*zgamma-ygamma*zeta) - &
-! ! xeta*(yxi*zgamma-ygamma*zxi) + &
-! ! xgamma*(yxi*zeta-yeta*zxi)
-
-! ! ! can ignore negative jacobian in mesher if needed when debugging code
-! ! if(jacobian <= ZERO) call exit_MPI(myrank,'3D Jacobian undefined')
-
-! ! ! invert the relation (Fletcher p. 50 vol. 2)
-! ! xix = (yeta*zgamma-ygamma*zeta) / jacobian
-! ! xiy = (xgamma*zeta-xeta*zgamma) / jacobian
-! ! xiz = (xeta*ygamma-xgamma*yeta) / jacobian
-! ! etax = (ygamma*zxi-yxi*zgamma) / jacobian
-! ! etay = (xxi*zgamma-xgamma*zxi) / jacobian
-! ! etaz = (xgamma*yxi-xxi*ygamma) / jacobian
-! ! gammax = (yxi*zeta-yeta*zxi) / jacobian
-! ! gammay = (xeta*zxi-xxi*zeta) / jacobian
-! ! gammaz = (xxi*yeta-xeta*yxi) / jacobian
-
-! ! ! compute and store the jacobian for the solver
-! ! jacobian = 1. / (xix*(etay*gammaz-etaz*gammay) &
-! ! -xiy*(etax*gammaz-etaz*gammax) &
-! ! +xiz*(etax*gammay-etay*gammax))
-
-! ! ! save the derivatives and the jacobian
-
-! ! ! distinguish between single and double precision for reals
-! ! if(CUSTOM_REAL == SIZE_REAL) then
-! ! xixstore(i,j,k,ispec) = sngl(xix)
-! ! xiystore(i,j,k,ispec) = sngl(xiy)
-! ! xizstore(i,j,k,ispec) = sngl(xiz)
-! ! etaxstore(i,j,k,ispec) = sngl(etax)
-! ! etaystore(i,j,k,ispec) = sngl(etay)
-! ! etazstore(i,j,k,ispec) = sngl(etaz)
-! ! gammaxstore(i,j,k,ispec) = sngl(gammax)
-! ! gammaystore(i,j,k,ispec) = sngl(gammay)
-! ! gammazstore(i,j,k,ispec) = sngl(gammaz)
-! ! jacobianstore(i,j,k,ispec) = sngl(jacobian)
-! ! else
-! ! xixstore(i,j,k,ispec) = xix
-! ! xiystore(i,j,k,ispec) = xiy
-! ! xizstore(i,j,k,ispec) = xiz
-! ! etaxstore(i,j,k,ispec) = etax
-! ! etaystore(i,j,k,ispec) = etay
-! ! etazstore(i,j,k,ispec) = etaz
-! ! gammaxstore(i,j,k,ispec) = gammax
-! ! gammaystore(i,j,k,ispec) = gammay
-! ! gammazstore(i,j,k,ispec) = gammaz
-! ! jacobianstore(i,j,k,ispec) = jacobian
-! ! endif
-
-! xstore(i,j,k,ispec) = xmesh
-! ystore(i,j,k,ispec) = ymesh
-! zstore(i,j,k,ispec) = zmesh
-
-! enddo
-! enddo
-! enddo
-
xstore(1,1,1,ispec) = xelm(1)
ystore(1,1,1,ispec) = yelm(1)
zstore(1,1,1,ispec) = zelm(1)
Modified: seismo/3D/SPECFEM3D/trunk/src/shared/create_movie_shakemap_AVS_DX_GMT.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/shared/create_movie_shakemap_AVS_DX_GMT.f90 2012-08-17 00:10:01 UTC (rev 20586)
+++ seismo/3D/SPECFEM3D/trunk/src/shared/create_movie_shakemap_AVS_DX_GMT.f90 2012-08-17 00:28:33 UTC (rev 20587)
@@ -409,9 +409,7 @@
do i = 1,NGLLX-1
ieoff = NGNOD2D_AVS_DX*(ielm+(i-1)+(j-1)*(NGLLX-1))
do ilocnum = 1,NGNOD2D_AVS_DX
- ! do k = 1,NGNOD2D_AVS_DX
-
if(ilocnum == 1) then
xp(ieoff+ilocnum) = dble(x(i,j))
yp(ieoff+ilocnum) = dble(y(i,j))
@@ -425,11 +423,6 @@
zp(ieoff+ilocnum) = dble(z(i+1,j+1))
field_display(ieoff+ilocnum) = dble(display(i+1,j+1))
- ! xp(ieoff+ilocnum) = dble(x(i+1,j))
- ! yp(ieoff+ilocnum) = dble(y(i+1,j))
- ! zp(ieoff+ilocnum) = dble(z(i+1,j))
- ! field_display(ieoff+ilocnum) = dble(display(i+1,j))
-
elseif(ilocnum == 3) then
! accounts for different ordering of square points
@@ -438,10 +431,6 @@
zp(ieoff+ilocnum) = dble(z(i+1,j))
field_display(ieoff+ilocnum) = dble(display(i+1,j))
- ! xp(ieoff+ilocnum) = dble(x(i+1,j+1))
- ! yp(ieoff+ilocnum) = dble(y(i+1,j+1))
- ! zp(ieoff+ilocnum) = dble(z(i+1,j+1))
- ! field_display(ieoff+ilocnum) = dble(display(i+1,j+1))
else
xp(ieoff+ilocnum) = dble(x(i,j+1))
yp(ieoff+ilocnum) = dble(y(i,j+1))
@@ -451,13 +440,6 @@
enddo
- !if( j==1 .and. ispec==1) then
- !print*,'p1',xp(ieoff+1),yp(ieoff+1),zp(ieoff+1)
- !print*,'p2',xp(ieoff+2),yp(ieoff+2),zp(ieoff+2)
- !print*,'p3',xp(ieoff+3),yp(ieoff+3),zp(ieoff+3)
- !print*,'p4',xp(ieoff+4),yp(ieoff+4),zp(ieoff+4)
- !endif
-
enddo
enddo
@@ -522,8 +504,6 @@
! velocity z-component
field_display(ilocnum+ieoff) = vectorz
endif
- ! takes norm of velocity vector
- !field_display(ilocnum+ieoff) =sqrt(vectorz**2+vectory**2+vectorx**2)
endif
enddo
@@ -706,8 +686,6 @@
if(USE_OPENDX) then
write(11,*) sngl(xp_save(ilocnum+ieoff)),sngl(yp_save(ilocnum+ieoff)),sngl(zp_save(ilocnum+ieoff))
else if(USE_AVS) then
- !write(11,'(i9,3f16.6)') ireorder(ibool_number),xp_save(ilocnum+ieoff), &
- ! yp_save(ilocnum+ieoff),zp_save(ilocnum+ieoff)
write(11,*) ireorder(ibool_number),sngl(xp_save(ilocnum+ieoff)), &
sngl(yp_save(ilocnum+ieoff)),sngl(zp_save(ilocnum+ieoff))
endif
@@ -759,14 +737,12 @@
if(plot_shaking_map) then
write(11,*) sngl(field_display(ilocnum+ieoff))
else
- !write(11,"(f7.2)") field_display(ilocnum+ieoff)
write(11,*) sngl(field_display(ilocnum+ieoff))
endif
else
if(plot_shaking_map) then
write(11,*) ireorder(ibool_number),sngl(field_display(ilocnum+ieoff))
else
- !write(11,"(i10,1x,f7.2)") ireorder(ibool_number),field_display(ilocnum+ieoff)
write(11,*) ireorder(ibool_number),sngl(field_display(ilocnum+ieoff))
endif
endif
@@ -802,7 +778,6 @@
if(USE_GMT) print *,'GMT files are stored in ', trim(OUTPUT_FILES), '/gmt_*.xyz'
print *
-
deallocate(store_val_x)
deallocate(store_val_y)
deallocate(store_val_z)
Modified: seismo/3D/SPECFEM3D/trunk/src/shared/get_element_face.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/shared/get_element_face.f90 2012-08-17 00:10:01 UTC (rev 20586)
+++ seismo/3D/SPECFEM3D/trunk/src/shared/get_element_face.f90 2012-08-17 00:28:33 UTC (rev 20587)
@@ -96,17 +96,6 @@
enddo
midpoint(:) = midpoint(:) / 4.0
- ! checks: this holds only for planar face
- !if( midpoint(1) /= (xcoord(1)+xcoord(3))/2.0 .or. midpoint(1) /= (xcoord(2)+xcoord(4))/2.0 ) then
- ! print*,'error midpoint x:',midpoint(1),(xcoord(1)+xcoord(3))/2.0,(xcoord(2)+xcoord(4))/2.0
- !endif
- !if( midpoint(2) /= (ycoord(1)+ycoord(3))/2.0 .or. midpoint(2) /= (ycoord(2)+ycoord(4))/2.0 ) then
- ! print*,'error midpoint y:',midpoint(1),(ycoord(1)+ycoord(3))/2.0,(ycoord(2)+ycoord(4))/2.0
- !endif
- !if( midpoint(3) /= (zcoord(1)+zcoord(3))/2.0 .or. midpoint(3) /= (zcoord(2)+zcoord(4))/2.0 ) then
- ! print*,'error midpoint z:',midpoint(1),(zcoord(1)+zcoord(3))/2.0,(zcoord(2)+zcoord(4))/2.0
- !endif
-
! determines element face by minimum distance of midpoints
midpoint_faces(:,:) = 0.0
do ifa=1,6
@@ -140,7 +129,6 @@
iloc = minloc(midpoint_distances)
! checks that found midpoint is close enough
- !print*,'face:', midpoint_distances(iloc(1))
if( midpoint_distances(iloc(1)) > 1.e-4 * &
( (xcoord(1)-xcoord(2))**2 &
+ (ycoord(1)-ycoord(2))**2 &
@@ -163,16 +151,6 @@
else
iface_id = iloc(1)
- !print*,'face:',iface_id
- !do icorner=1,NGNOD2D
- ! 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))
- ! iglob = ibool(i,j,k,ispec)
- ! print*,'corner:',icorner,'xyz:',sngl(xstore_dummy(iglob)), &
- ! sngl(ystore_dummy(iglob)),sngl(zstore_dummy(iglob))
- !enddo
-
endif
end subroutine get_element_face_id
@@ -195,7 +173,6 @@
integer :: NGLLA,NGLLB
integer,dimension(3,NGLLA,NGLLB) :: ijk_face
-! integer :: icorner,i,j,k,iglob,iloc(1)
integer :: i,j,k
integer :: ngll,i_gll,j_gll,k_gll
@@ -291,80 +268,6 @@
print*,'error element face ngll:',ngll,NGLLA,NGLLB
stop 'error element face ngll'
endif
-!
-!! corner locations
-! do icorner=1,NGNOD2D
-! i = iface_all_corner_ijk(1,icorner,iface)
-! j = iface_all_corner_ijk(2,icorner,iface)
-! k = iface_all_corner_ijk(3,icorner,iface)
-! iglob = ibool(i,j,k,ispec)
-! xcoord_iboun(icorner) = xstore_dummy(iglob)
-! ycoord_iboun(icorner) = ystore_dummy(iglob)
-! zcoord_iboun(icorner) = zstore_dummy(iglob)
-! ! looks at values
-! !print*,'corner:',icorner,'xyz:',sngl(xcoord_iboun(icorner)),sngl(ycoord_iboun(icorner)),sngl(zcoord_iboun(icorner))
-! enddo
-!
-!! determines initial orientation given by three corners of the face
-! ! (CUBIT orders corners such that normal points outwards of element)
-! ! cross-product of vectors from corner 1 to corner 2 and from corner 1 to corner 3
-! face_n(1) = (ycoord(2)-ycoord(1))*(zcoord(3)-zcoord(1)) - (zcoord(2)-zcoord(1))*(ycoord(3)-ycoord(1))
-! face_n(2) = - (xcoord(2)-xcoord(1))*(zcoord(3)-zcoord(1)) + (zcoord(2)-zcoord(1))*(xcoord(3)-xcoord(1))
-! face_n(3) = (xcoord(2)-xcoord(1))*(ycoord(3)-ycoord(1)) - (ycoord(2)-ycoord(1))*(xcoord(3)-xcoord(1))
-! face_n(:) = face_n(:)/(sqrt( face_n(1)**2 + face_n(2)**2 + face_n(3)**2) )
-!
-!! checks that this normal direction is outwards of element:
-! ! takes additional corner out of face plane and determines scalarproduct to normal
-! select case( iface )
-! case(1) ! opposite to xmin face
-! iglob = ibool(NGLLX,1,1,ispec)
-! case(2) ! opposite to xmax face
-! iglob = ibool(1,1,1,ispec)
-! case(3) ! opposite to ymin face
-! iglob = ibool(1,NGLLY,1,ispec)
-! case(4) ! opposite to ymax face
-! iglob = ibool(1,1,1,ispec)
-! case(5) ! opposite to bottom
-! iglob = ibool(1,1,NGLLZ,ispec)
-! case(6) ! opposite to top
-! iglob = ibool(1,1,1,ispec)
-! end select
-! ! vector from corner 1 to this opposite one
-! xcoord(4) = xstore_dummy(iglob) - xcoord(1)
-! ycoord(4) = ystore_dummy(iglob) - ycoord(1)
-! zcoord(4) = zstore_dummy(iglob) - zcoord(1)
-!
-! ! scalarproduct
-! tmp = xcoord(4)*face_n(1) + ycoord(4)*face_n(2) + zcoord(4)*face_n(3)
-!
-! ! makes sure normal points outwards, that is points away from this additional corner and scalarproduct is negative
-! if( tmp > 0.0 ) then
-! face_n(:) = - face_n(:)
-! endif
-! !print*,'face ',iface,'scalarproduct:',tmp
-!
-!! determines orientation of gll corner locations and sets it such that normal points outwards
-! ! cross-product
-! face_ntmp(1) = (ycoord_iboun(2)-ycoord_iboun(1))*(zcoord_iboun(3)-zcoord_iboun(1)) &
-! - (zcoord_iboun(2)-zcoord_iboun(1))*(ycoord_iboun(3)-ycoord_iboun(1))
-! face_ntmp(2) = - (xcoord_iboun(2)-xcoord_iboun(1))*(zcoord_iboun(3)-zcoord_iboun(1)) &
-! + (zcoord_iboun(2)-zcoord_iboun(1))*(xcoord_iboun(3)-xcoord_iboun(1))
-! face_ntmp(3) = (xcoord_iboun(2)-xcoord_iboun(1))*(ycoord_iboun(3)-ycoord_iboun(1))&
-! - (ycoord_iboun(2)-ycoord_iboun(1))*(xcoord_iboun(3)-xcoord_iboun(1))
-! face_ntmp(:) = face_ntmp(:)/(sqrt( face_ntmp(1)**2 + face_ntmp(2)**2 + face_ntmp(3)**2) )
-! if( abs( (face_n(1)-face_ntmp(1))**2+(face_n(2)-face_ntmp(2))**2+(face_n(3)-face_ntmp(3))**2) > 0.1 ) then
-! !print*,'error orientation face 1:',ispec,face_n(:)
-! !swap corners 2 and 4 ( switches clockwise / anti-clockwise )
-! tmp = xcoord_iboun(2)
-! xcoord_iboun(2) = xcoord_iboun(4)
-! xcoord_iboun(4) = tmp
-! tmp = ycoord_iboun(2)
-! ycoord_iboun(2) = ycoord_iboun(4)
-! ycoord_iboun(4) = tmp
-! tmp = zcoord_iboun(2)
-! zcoord_iboun(2) = zcoord_iboun(4)
-! zcoord_iboun(4) = tmp
-! endif
end subroutine get_element_face_gll_indices
@@ -401,7 +304,6 @@
real(kind=CUSTOM_REAL) :: face_n(3),tmp,v_tmp(3)
integer :: iglob
-
! determines initial orientation given by three corners on the face
! cross-product of vectors from corner 1 to corner 2 and from corner 1 to corner 3
face_n(1) = (ycoord(2)-ycoord(1))*(zcoord(3)-zcoord(1)) - (zcoord(2)-zcoord(1))*(ycoord(3)-ycoord(1))
@@ -456,13 +358,9 @@
! 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 ) then
- !print*,'element face normal: orientation ',ispec,iface,tmp
- !print*,'face normal: ',face_n(:)
- !print*,' normal: ',normal(:)
!swap
normal(:) = - normal(:)
endif
- !print*,'face ',iface,'scalarproduct:',tmp
end subroutine get_element_face_normal
@@ -613,3 +511,4 @@
enddo
end subroutine get_element_corners
+
Modified: seismo/3D/SPECFEM3D/trunk/src/shared/get_jacobian_boundaries.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/shared/get_jacobian_boundaries.f90 2012-08-17 00:10:01 UTC (rev 20586)
+++ seismo/3D/SPECFEM3D/trunk/src/shared/get_jacobian_boundaries.f90 2012-08-17 00:28:33 UTC (rev 20587)
@@ -24,7 +24,6 @@
!
!=====================================================================
-
subroutine get_jacobian_boundary_face(myrank,nspec, &
xstore_dummy,ystore_dummy,zstore_dummy,ibool,nglob,&
dershape2D_x,dershape2D_y,dershape2D_bottom,dershape2D_top, &
@@ -262,7 +261,7 @@
! This subroutine recompute the 3D jacobian for one element
! based upon 125 GLL points
-! Hejun Zhu OCT16,2009
+! Hejun Zhu, Oct 16, 2009
! input: myrank,
! xstore,ystore,zstore ----- input position
@@ -407,527 +406,3 @@
end subroutine recalc_jacobian_gll2D
-!
-!------------------------------------------------------------------------------------------------
-!
-!
-! subroutine get_jacobian_boundaries(myrank,iboun,nspec, &
-! xstore_dummy,ystore_dummy,zstore_dummy,ibool,nglob,&
-! dershape2D_x,dershape2D_y,dershape2D_bottom,dershape2D_top, &
-! wgllwgll_xy,wgllwgll_xz,wgllwgll_yz,&
-! ibelm_xmin,ibelm_xmax,ibelm_ymin,ibelm_ymax,ibelm_bottom,ibelm_top, &
-! xcoord_iboun,ycoord_iboun,zcoord_iboun, &
-! nspec2D_xmin,nspec2D_xmax,nspec2D_ymin,nspec2D_ymax, &
-! jacobian2D_xmin,jacobian2D_xmax, &
-! jacobian2D_ymin,jacobian2D_ymax, &
-! jacobian2D_bottom,jacobian2D_top, &
-! normal_xmin,normal_xmax, &
-! normal_ymin,normal_ymax, &
-! normal_bottom,normal_top, &
-! NSPEC2D_BOTTOM,NSPEC2D_TOP)
-!
-! implicit none
-!
-! include "constants.h"
-!
-! integer nspec,myrank,nglob
-!
-!! arrays with the mesh
-! integer, dimension(NGLLX,NGLLY,NGLLZ,nspec) :: ibool
-! real(kind=CUSTOM_REAL) :: xstore_dummy(nglob),ystore_dummy(nglob),zstore_dummy(nglob)
-!
-!
-!! absorbing boundaries
-!! (careful with array bounds, no need for NSPEC2DMAX_XMIN_XMAX & NSPEC2DMAX_YMIN_YMAX anymore)
-! integer :: nspec2D_xmin, nspec2D_xmax, nspec2D_ymin, nspec2D_ymax, NSPEC2D_BOTTOM, NSPEC2D_TOP
-! integer, dimension(nspec2D_xmin) :: ibelm_xmin
-! integer, dimension(nspec2D_xmax) :: ibelm_xmax
-! integer, dimension(nspec2D_ymin) :: ibelm_ymin
-! integer, dimension(nspec2D_ymax) :: ibelm_ymax
-! integer, dimension(NSPEC2D_BOTTOM) :: ibelm_bottom
-! integer, dimension(NSPEC2D_TOP) :: ibelm_top
-!
-! logical iboun(6,nspec)
-! real(kind=CUSTOM_REAL), dimension(NGNOD2D,6,nspec) :: xcoord_iboun,ycoord_iboun,zcoord_iboun
-!
-!! double precision xstore(NGLLX,NGLLY,NGLLZ,nspec)
-!! double precision ystore(NGLLX,NGLLY,NGLLZ,nspec)
-!! double precision zstore(NGLLX,NGLLY,NGLLZ,nspec)
-!
-! real(kind=CUSTOM_REAL) jacobian2D_xmin(NGLLY,NGLLZ,NSPEC2D_xmin)
-! real(kind=CUSTOM_REAL) jacobian2D_xmax(NGLLY,NGLLZ,NSPEC2D_xmax)
-! real(kind=CUSTOM_REAL) jacobian2D_ymin(NGLLX,NGLLZ,NSPEC2D_ymin)
-! real(kind=CUSTOM_REAL) jacobian2D_ymax(NGLLX,NGLLZ,NSPEC2D_ymax)
-! real(kind=CUSTOM_REAL) jacobian2D_bottom(NGLLX,NGLLY,NSPEC2D_BOTTOM)
-! real(kind=CUSTOM_REAL) jacobian2D_top(NGLLX,NGLLY,NSPEC2D_TOP)
-!
-! real(kind=CUSTOM_REAL) normal_xmin(NDIM,NGLLY,NGLLZ,NSPEC2D_xmin)
-! real(kind=CUSTOM_REAL) normal_xmax(NDIM,NGLLY,NGLLZ,NSPEC2D_xmax)
-! real(kind=CUSTOM_REAL) normal_ymin(NDIM,NGLLX,NGLLZ,NSPEC2D_ymin)
-! real(kind=CUSTOM_REAL) normal_ymax(NDIM,NGLLX,NGLLZ,NSPEC2D_ymax)
-! real(kind=CUSTOM_REAL) normal_bottom(NDIM,NGLLX,NGLLY,NSPEC2D_BOTTOM)
-! real(kind=CUSTOM_REAL) normal_top(NDIM,NGLLX,NGLLY,NSPEC2D_TOP)
-!
-! double precision dershape2D_x(NDIM2D,NGNOD2D,NGLLY,NGLLZ)
-! double precision dershape2D_y(NDIM2D,NGNOD2D,NGLLX,NGLLZ)
-! double precision dershape2D_bottom(NDIM2D,NGNOD2D,NGLLX,NGLLY)
-! double precision dershape2D_top(NDIM2D,NGNOD2D,NGLLX,NGLLY)
-!
-! double precision, dimension(NGLLX,NGLLY) :: wgllwgll_xy
-! double precision, dimension(NGLLX,NGLLZ) :: wgllwgll_xz
-! double precision, dimension(NGLLY,NGLLZ) :: wgllwgll_yz
-!
-! double precision xelm(NGNOD2D),yelm(NGNOD2D),zelm(NGNOD2D)
-!
-!! element numbering
-! integer ispec,i,j
-!
-!! counters to keep track of number of elements on each of the boundaries
-! integer ispecb1,ispecb2,ispecb3,ispecb4,ispecb5,ispecb6
-!
-!
-!! check that the parameter file is correct
-! if(NGNOD /= 8) call exit_MPI(myrank,'elements should have 8 control nodes')
-! if(NGNOD2D /= 4) call exit_MPI(myrank,'surface elements should have 4 control nodes')
-!
-! ispecb1 = 0
-! ispecb2 = 0
-! ispecb3 = 0
-! ispecb4 = 0
-! ispecb5 = 0
-! ispecb6 = 0
-!
-! do ispec=1,nspec
-!
-!! determine if the element falls on a boundary
-!
-!! on boundary: xmin
-!
-! if(iboun(1,ispec)) then
-!
-! ispecb1=ispecb1+1
-! ibelm_xmin(ispecb1)=ispec
-!
-!! specify the 4 nodes for the 2-D boundary element
-!! i.e. face (0,0,0),(0,1,0),(0,1,1),(0,0,1)
-!
-!! careful: these points may not be on the xmin face for unstructured grids
-!! xelm(1)=xstore(1,1,1,ispec)
-!! yelm(1)=ystore(1,1,1,ispec)
-!! zelm(1)=zstore(1,1,1,ispec)
-!! xelm(2)=xstore(1,NGLLY,1,ispec)
-!! yelm(2)=ystore(1,NGLLY,1,ispec)
-!! zelm(2)=zstore(1,NGLLY,1,ispec)
-!! xelm(3)=xstore(1,NGLLY,NGLLZ,ispec)
-!! yelm(3)=ystore(1,NGLLY,NGLLZ,ispec)
-!! zelm(3)=zstore(1,NGLLY,NGLLZ,ispec)
-!! xelm(4)=xstore(1,1,NGLLZ,ispec)
-!! yelm(4)=ystore(1,1,NGLLZ,ispec)
-!! zelm(4)=zstore(1,1,NGLLZ,ispec)
-!
-! xelm(1)=xstore_dummy( ibool(1,1,1,ispec) )
-! yelm(1)=ystore_dummy( ibool(1,1,1,ispec) )
-! zelm(1)=zstore_dummy( ibool(1,1,1,ispec) )
-! xelm(2)=xstore_dummy( ibool(1,NGLLY,1,ispec) )
-! yelm(2)=ystore_dummy( ibool(1,NGLLY,1,ispec) )
-! zelm(2)=zstore_dummy( ibool(1,NGLLY,1,ispec) )
-! xelm(3)=xstore_dummy( ibool(1,NGLLY,NGLLZ,ispec) )
-! yelm(3)=ystore_dummy( ibool(1,NGLLY,NGLLZ,ispec) )
-! zelm(3)=zstore_dummy( ibool(1,NGLLY,NGLLZ,ispec) )
-! xelm(4)=xstore_dummy( ibool(1,1,NGLLZ,ispec) )
-! yelm(4)=ystore_dummy( ibool(1,1,NGLLZ,ispec) )
-! zelm(4)=zstore_dummy( ibool(1,1,NGLLZ,ispec) )
-!
-!! takes coordinates from boundary faces
-!! do i=1,NGNOD2D
-!! xelm(i) = xcoord_iboun(i,1,ispec)
-!! yelm(i) = ycoord_iboun(i,1,ispec)
-!! zelm(i) = zcoord_iboun(i,1,ispec)
-!! enddo
-!
-! call compute_jacobian_2D(myrank,ispecb1,xelm,yelm,zelm, &
-! dershape2D_x,wgllwgll_yz, &
-! jacobian2D_xmin,normal_xmin,NGLLY,NGLLZ,NSPEC2D_xmin)
-!
-! ! normal convention: points away from element
-! ! switches normal direction if necessary
-! do i=1,NGLLY
-! do j=1,NGLLZ
-! call get_element_face_normal(ispecb1, 1, xelm,yelm,zelm, &
-! ibool,nspec,nglob, &
-! xstore_dummy,ystore_dummy,zstore_dummy, &
-! normal_xmin(:,i,j,ispecb1) )
-! enddo
-! enddo
-!
-! endif
-!
-!! on boundary: xmax
-!
-! if(iboun(2,ispec)) then
-!
-! ispecb2=ispecb2+1
-! ibelm_xmax(ispecb2)=ispec
-!
-!! careful...
-!! specify the 4 nodes for the 2-D boundary element
-!! xelm(1)=xstore(NGLLX,1,1,ispec)
-!! yelm(1)=ystore(NGLLX,1,1,ispec)
-!! zelm(1)=zstore(NGLLX,1,1,ispec)
-!! xelm(2)=xstore(NGLLX,NGLLY,1,ispec)
-!! yelm(2)=ystore(NGLLX,NGLLY,1,ispec)
-!! zelm(2)=zstore(NGLLX,NGLLY,1,ispec)
-!! xelm(3)=xstore(NGLLX,NGLLY,NGLLZ,ispec)
-!! yelm(3)=ystore(NGLLX,NGLLY,NGLLZ,ispec)
-!! zelm(3)=zstore(NGLLX,NGLLY,NGLLZ,ispec)
-!! xelm(4)=xstore(NGLLX,1,NGLLZ,ispec)
-!! yelm(4)=ystore(NGLLX,1,NGLLZ,ispec)
-!! zelm(4)=zstore(NGLLX,1,NGLLZ,ispec)
-!
-! xelm(1)=xstore_dummy( ibool(NGLLX,1,1,ispec) )
-! yelm(1)=ystore_dummy( ibool(NGLLX,1,1,ispec) )
-! zelm(1)=zstore_dummy( ibool(NGLLX,1,1,ispec) )
-! xelm(2)=xstore_dummy( ibool(NGLLX,NGLLY,1,ispec) )
-! yelm(2)=ystore_dummy( ibool(NGLLX,NGLLY,1,ispec) )
-! zelm(2)=zstore_dummy( ibool(NGLLX,NGLLY,1,ispec) )
-! xelm(3)=xstore_dummy( ibool(NGLLX,NGLLY,NGLLZ,ispec) )
-! yelm(3)=ystore_dummy( ibool(NGLLX,NGLLY,NGLLZ,ispec) )
-! zelm(3)=zstore_dummy( ibool(NGLLX,NGLLY,NGLLZ,ispec) )
-! xelm(4)=xstore_dummy( ibool(NGLLX,1,NGLLZ,ispec) )
-! yelm(4)=ystore_dummy( ibool(NGLLX,1,NGLLZ,ispec) )
-! zelm(4)=zstore_dummy( ibool(NGLLX,1,NGLLZ,ispec) )
-!
-!! takes coordinates from boundary faces
-!! do i=1,NGNOD2D
-!! xelm(i) = xcoord_iboun(i,2,ispec)
-!! yelm(i) = ycoord_iboun(i,2,ispec)
-!! zelm(i) = zcoord_iboun(i,2,ispec)
-!! enddo
-!
-! call compute_jacobian_2D(myrank,ispecb2,xelm,yelm,zelm, &
-! dershape2D_x,wgllwgll_yz, &
-! jacobian2D_xmax,normal_xmax,NGLLY,NGLLZ,NSPEC2D_xmax)
-!
-! ! normal convention: points away from element
-! ! switch normal direction if necessary
-! do i=1,NGLLY
-! do j=1,NGLLZ
-! call get_element_face_normal(ispecb2, 2, xelm,yelm,zelm, &
-! ibool,nspec,nglob, &
-! xstore_dummy,ystore_dummy,zstore_dummy, &
-! normal_xmax(:,i,j,ispecb2) )
-! enddo
-! enddo
-!
-! endif
-!
-!! on boundary: ymin
-!
-! if(iboun(3,ispec)) then
-!
-! ispecb3=ispecb3+1
-! ibelm_ymin(ispecb3)=ispec
-!
-!! careful...
-!! specify the 4 nodes for the 2-D boundary element
-!! xelm(1)=xstore(1,1,1,ispec)
-!! yelm(1)=ystore(1,1,1,ispec)
-!! zelm(1)=zstore(1,1,1,ispec)
-!! xelm(2)=xstore(NGLLX,1,1,ispec)
-!! yelm(2)=ystore(NGLLX,1,1,ispec)
-!! zelm(2)=zstore(NGLLX,1,1,ispec)
-!! xelm(3)=xstore(NGLLX,1,NGLLZ,ispec)
-!! yelm(3)=ystore(NGLLX,1,NGLLZ,ispec)
-!! zelm(3)=zstore(NGLLX,1,NGLLZ,ispec)
-!! xelm(4)=xstore(1,1,NGLLZ,ispec)
-!! yelm(4)=ystore(1,1,NGLLZ,ispec)
-!! zelm(4)=zstore(1,1,NGLLZ,ispec)
-!
-! xelm(1)=xstore_dummy( ibool(1,1,1,ispec) )
-! yelm(1)=ystore_dummy( ibool(1,1,1,ispec) )
-! zelm(1)=zstore_dummy( ibool(1,1,1,ispec) )
-! xelm(2)=xstore_dummy( ibool(NGLLX,1,1,ispec) )
-! yelm(2)=ystore_dummy( ibool(NGLLX,1,1,ispec) )
-! zelm(2)=zstore_dummy( ibool(NGLLX,1,1,ispec) )
-! xelm(3)=xstore_dummy( ibool(NGLLX,1,NGLLZ,ispec) )
-! yelm(3)=ystore_dummy( ibool(NGLLX,1,NGLLZ,ispec) )
-! zelm(3)=zstore_dummy( ibool(NGLLX,1,NGLLZ,ispec) )
-! xelm(4)=xstore_dummy( ibool(1,1,NGLLZ,ispec) )
-! yelm(4)=ystore_dummy( ibool(1,1,NGLLZ,ispec) )
-! zelm(4)=zstore_dummy( ibool(1,1,NGLLZ,ispec) )
-!
-!! takes coordinates from boundary faces
-!! do i=1,NGNOD2D
-!! xelm(i) = xcoord_iboun(i,3,ispec)
-!! yelm(i) = ycoord_iboun(i,3,ispec)
-!! zelm(i) = zcoord_iboun(i,3,ispec)
-!! enddo
-!
-! call compute_jacobian_2D(myrank,ispecb3,xelm,yelm,zelm, &
-! dershape2D_y,wgllwgll_xz, &
-! jacobian2D_ymin,normal_ymin,NGLLX,NGLLZ,NSPEC2D_ymin)
-!
-! ! normal convention: points away from element
-! ! switch normal direction if necessary
-! do i=1,NGLLX
-! do j=1,NGLLZ
-! call get_element_face_normal(ispecb3, 3, xelm,yelm,zelm, &
-! ibool,nspec,nglob, &
-! xstore_dummy,ystore_dummy,zstore_dummy, &
-! normal_ymin(:,i,j,ispecb3) )
-! enddo
-! enddo
-!
-!
-! endif
-!
-!! on boundary: ymax
-!
-! if(iboun(4,ispec)) then
-!
-! ispecb4=ispecb4+1
-! ibelm_ymax(ispecb4)=ispec
-!
-!!careful...
-!! specify the 4 nodes for the 2-D boundary element
-!! xelm(1)=xstore(1,NGLLY,1,ispec)
-!! yelm(1)=ystore(1,NGLLY,1,ispec)
-!! zelm(1)=zstore(1,NGLLY,1,ispec)
-!! xelm(2)=xstore(NGLLX,NGLLY,1,ispec)
-!! yelm(2)=ystore(NGLLX,NGLLY,1,ispec)
-!! zelm(2)=zstore(NGLLX,NGLLY,1,ispec)
-!! xelm(3)=xstore(NGLLX,NGLLY,NGLLZ,ispec)
-!! yelm(3)=ystore(NGLLX,NGLLY,NGLLZ,ispec)
-!! zelm(3)=zstore(NGLLX,NGLLY,NGLLZ,ispec)
-!! xelm(4)=xstore(1,NGLLY,NGLLZ,ispec)
-!! yelm(4)=ystore(1,NGLLY,NGLLZ,ispec)
-!! zelm(4)=zstore(1,NGLLY,NGLLZ,ispec)
-!
-! xelm(1)=xstore_dummy( ibool(1,NGLLY,1,ispec) )
-! yelm(1)=ystore_dummy( ibool(1,NGLLY,1,ispec) )
-! zelm(1)=zstore_dummy( ibool(1,NGLLY,1,ispec) )
-! xelm(2)=xstore_dummy( ibool(NGLLX,NGLLY,1,ispec) )
-! yelm(2)=ystore_dummy( ibool(NGLLX,NGLLY,1,ispec) )
-! zelm(2)=zstore_dummy( ibool(NGLLX,NGLLY,1,ispec) )
-! xelm(3)=xstore_dummy( ibool(NGLLX,NGLLY,NGLLZ,ispec) )
-! yelm(3)=ystore_dummy( ibool(NGLLX,NGLLY,NGLLZ,ispec) )
-! zelm(3)=zstore_dummy( ibool(NGLLX,NGLLY,NGLLZ,ispec) )
-! xelm(4)=xstore_dummy( ibool(1,NGLLY,NGLLZ,ispec) )
-! yelm(4)=ystore_dummy( ibool(1,NGLLY,NGLLZ,ispec) )
-! zelm(4)=zstore_dummy( ibool(1,NGLLY,NGLLZ,ispec) )
-!
-!! takes coordinates from boundary faces
-!! do i=1,NGNOD2D
-!! xelm(i) = xcoord_iboun(i,4,ispec)
-!! yelm(i) = ycoord_iboun(i,4,ispec)
-!! zelm(i) = zcoord_iboun(i,4,ispec)
-!! enddo
-!!
-! call compute_jacobian_2D(myrank,ispecb4,xelm,yelm,zelm, &
-! dershape2D_y, wgllwgll_xz, &
-! jacobian2D_ymax,normal_ymax,NGLLX,NGLLZ,NSPEC2D_ymax)
-!
-! ! normal convention: points away from element
-! ! switch normal direction if necessary
-! do i=1,NGLLX
-! do j=1,NGLLZ
-! call get_element_face_normal(ispecb4, 4, xelm,yelm,zelm, &
-! ibool,nspec,nglob, &
-! xstore_dummy,ystore_dummy,zstore_dummy, &
-! normal_ymax(:,i,j,ispecb4) )
-! enddo
-! enddo
-!
-! endif
-!
-!! on boundary: bottom
-!
-! if(iboun(5,ispec)) then
-!
-! ispecb5=ispecb5+1
-! ibelm_bottom(ispecb5)=ispec
-!
-!! careful...
-!! for bottom, this might be actually working... when mesh is oriented along z direction...
-!! xelm(1)=xstore(1,1,1,ispec)
-!! yelm(1)=ystore(1,1,1,ispec)
-!! zelm(1)=zstore(1,1,1,ispec)
-!! xelm(2)=xstore(NGLLX,1,1,ispec)
-!! yelm(2)=ystore(NGLLX,1,1,ispec)
-!! zelm(2)=zstore(NGLLX,1,1,ispec)
-!! xelm(3)=xstore(NGLLX,NGLLY,1,ispec)
-!! yelm(3)=ystore(NGLLX,NGLLY,1,ispec)
-!! zelm(3)=zstore(NGLLX,NGLLY,1,ispec)
-!! xelm(4)=xstore(1,NGLLY,1,ispec)
-!! yelm(4)=ystore(1,NGLLY,1,ispec)
-!! zelm(4)=zstore(1,NGLLY,1,ispec)
-!
-! xelm(1)=xstore_dummy( ibool(1,1,1,ispec) )
-! yelm(1)=ystore_dummy( ibool(1,1,1,ispec) )
-! zelm(1)=zstore_dummy( ibool(1,1,1,ispec) )
-! xelm(2)=xstore_dummy( ibool(NGLLX,1,1,ispec) )
-! yelm(2)=ystore_dummy( ibool(NGLLX,1,1,ispec) )
-! zelm(2)=zstore_dummy( ibool(NGLLX,1,1,ispec) )
-! xelm(3)=xstore_dummy( ibool(NGLLX,NGLLY,1,ispec) )
-! yelm(3)=ystore_dummy( ibool(NGLLX,NGLLY,1,ispec) )
-! zelm(3)=zstore_dummy( ibool(NGLLX,NGLLY,1,ispec) )
-! xelm(4)=xstore_dummy( ibool(1,NGLLY,1,ispec) )
-! yelm(4)=ystore_dummy( ibool(1,NGLLY,1,ispec) )
-! zelm(4)=zstore_dummy( ibool(1,NGLLY,1,ispec) )
-!
-!
-!! takes coordinates from boundary faces
-!! do i=1,NGNOD2D
-!! xelm(i) = xcoord_iboun(i,5,ispec)
-!! yelm(i) = ycoord_iboun(i,5,ispec)
-!! zelm(i) = zcoord_iboun(i,5,ispec)
-!! enddo
-!
-! call compute_jacobian_2D(myrank,ispecb5,xelm,yelm,zelm,&
-! dershape2D_bottom,wgllwgll_xy, &
-! jacobian2D_bottom,normal_bottom,NGLLX,NGLLY,NSPEC2D_BOTTOM)
-!
-! ! normal convention: points away from element
-! ! switch normal direction if necessary
-! do i=1,NGLLX
-! do j=1,NGLLY
-! call get_element_face_normal(ispecb5, 5, xelm,yelm,zelm, &
-! ibool,nspec,nglob, &
-! xstore_dummy,ystore_dummy,zstore_dummy, &
-! normal_bottom(:,i,j,ispecb5) )
-! enddo
-! enddo
-!
-! endif
-!
-!! on boundary: top
-!
-! if(iboun(6,ispec)) then
-!
-! ispecb6=ispecb6+1
-! ibelm_top(ispecb6)=ispec
-!
-!! careful...
-!! for top, this might be working as well ... when mesh is oriented along z direction...
-!! xelm(1)=xstore(1,1,NGLLZ,ispec)
-!! yelm(1)=ystore(1,1,NGLLZ,ispec)
-!! zelm(1)=zstore(1,1,NGLLZ,ispec)
-!! xelm(2)=xstore(NGLLX,1,NGLLZ,ispec)
-!! yelm(2)=ystore(NGLLX,1,NGLLZ,ispec)
-!! zelm(2)=zstore(NGLLX,1,NGLLZ,ispec)
-!! xelm(3)=xstore(NGLLX,NGLLY,NGLLZ,ispec)
-!! yelm(3)=ystore(NGLLX,NGLLY,NGLLZ,ispec)
-!! zelm(3)=zstore(NGLLX,NGLLY,NGLLZ,ispec)
-!! xelm(4)=xstore(1,NGLLY,NGLLZ,ispec)
-!! yelm(4)=ystore(1,NGLLY,NGLLZ,ispec)
-!! zelm(4)=zstore(1,NGLLY,NGLLZ,ispec)
-!
-!
-!! takes coordinates from boundary faces
-!! do i=1,NGNOD2D
-!! xelm(i) = xcoord_iboun(i,6,ispec)
-!! yelm(i) = ycoord_iboun(i,6,ispec)
-!! zelm(i) = zcoord_iboun(i,6,ispec)
-!! enddo
-!
-! call compute_jacobian_2D(myrank,ispecb6,xelm,yelm,zelm,&
-! dershape2D_top, wgllwgll_xy, &
-! jacobian2D_top,normal_top,NGLLX,NGLLY,NSPEC2D_TOP)
-!
-! ! normal convention: points away from element
-! ! switch normal direction if necessary
-! do i=1,NGLLX
-! do j=1,NGLLY
-! call get_element_face_normal(ispecb6, 6, xelm,yelm,zelm, &
-! ibool,nspec,nglob, &
-! xstore_dummy,ystore_dummy,zstore_dummy, &
-! normal_top(:,i,j,ispecb6) )
-! enddo
-! enddo
-!
-! endif
-!
-! enddo
-!
-!! check theoretical value of elements
-!! if(ispecb1 /= NSPEC2D_xmin) call exit_MPI(myrank,'ispecb1 should equal NSPEC2D_xmin')
-!! if(ispecb2 /= NSPEC2D_xmax) call exit_MPI(myrank,'ispecb2 should equal NSPEC2D_xmax')
-!! if(ispecb3 /= NSPEC2D_ymin) call exit_MPI(myrank,'ispecb3 should equal NSPEC2D_ymin')
-!! if(ispecb4 /= NSPEC2D_ymax) call exit_MPI(myrank,'ispecb4 should equal NSPEC2D_ymax')
-!! if(ispecb5 /= NSPEC2D_BOTTOM) call exit_MPI(myrank,'ispecb5 should equal NSPEC2D_BOTTOM')
-!! if(ispecb6 /= NSPEC2D_TOP) call exit_MPI(myrank,'ispecb6 should equal NSPEC2D_TOP')
-!
-! end subroutine get_jacobian_boundaries
-!
-!! -------------------------------------------------------
-!
-! subroutine compute_jacobian_2D(myrank,ispecb,xelm,yelm,zelm, &
-! dershape2D,wgllwgll, &
-! jacobian2D,normal, &
-! NGLLA,NGLLB,NSPEC2DMAX_AB)
-!
-! implicit none
-!
-! include "constants.h"
-!
-!! generic routine that accepts any polynomial degree in each direction
-!
-! integer ispecb,NGLLA,NGLLB,NSPEC2DMAX_AB,myrank
-!
-! double precision xelm(NGNOD2D),yelm(NGNOD2D),zelm(NGNOD2D)
-! double precision dershape2D(NDIM2D,NGNOD2D,NGLLA,NGLLB)
-! double precision wgllwgll
-!
-! real(kind=CUSTOM_REAL) jacobian2D(NGLLA,NGLLB,NSPEC2DMAX_AB)
-! real(kind=CUSTOM_REAL) normal(3,NGLLA,NGLLB,NSPEC2DMAX_AB)
-!
-! integer i,j,ia
-! double precision xxi,xeta,yxi,yeta,zxi,zeta
-! double precision unx,uny,unz,jacobian
-!
-! do j=1,NGLLB
-! do i=1,NGLLA
-!
-! xxi=ZERO
-! xeta=ZERO
-! yxi=ZERO
-! yeta=ZERO
-! zxi=ZERO
-! zeta=ZERO
-! do ia=1,NGNOD2D
-! xxi=xxi+dershape2D(1,ia,i,j)*xelm(ia)
-! xeta=xeta+dershape2D(2,ia,i,j)*xelm(ia)
-! yxi=yxi+dershape2D(1,ia,i,j)*yelm(ia)
-! yeta=yeta+dershape2D(2,ia,i,j)*yelm(ia)
-! zxi=zxi+dershape2D(1,ia,i,j)*zelm(ia)
-! zeta=zeta+dershape2D(2,ia,i,j)*zelm(ia)
-! enddo
-!
-!! 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
-! jacobian2D(i,j,ispecb) = sngl(jacobian * wgllwgll(i,j) )
-! normal(1,i,j,ispecb)=sngl(unx/jacobian)
-! normal(2,i,j,ispecb)=sngl(uny/jacobian)
-! normal(3,i,j,ispecb)=sngl(unz/jacobian)
-! else
-! jacobian2D(i,j,ispecb) = jacobian * wgllwgll(i,j)
-! normal(1,i,j,ispecb)=unx/jacobian
-! normal(2,i,j,ispecb)=uny/jacobian
-! normal(3,i,j,ispecb)=unz/jacobian
-! endif
-!
-! enddo
-! enddo
-!
-! end subroutine compute_jacobian_2D
-!
-
More information about the CIG-COMMITS
mailing list