[cig-commits] r12935 - seismo/3D/SPECFEM3D/branches/update_temporary
nlegoff at geodynamics.org
nlegoff at geodynamics.org
Mon Sep 22 08:01:36 PDT 2008
Author: nlegoff
Date: 2008-09-22 08:01:36 -0700 (Mon, 22 Sep 2008)
New Revision: 12935
Added:
seismo/3D/SPECFEM3D/branches/update_temporary/sort_array_coordinates.f90
Modified:
seismo/3D/SPECFEM3D/branches/update_temporary/Makefile.in
seismo/3D/SPECFEM3D/branches/update_temporary/constants.h.in
seismo/3D/SPECFEM3D/branches/update_temporary/create_regions_mesh.f90
seismo/3D/SPECFEM3D/branches/update_temporary/meshfem3D.f90
Log:
added external mesh to meshfem3D.
Still lacks the generation of database (the partitioning), sources/receivers management and Stacey conditions.
Modified: seismo/3D/SPECFEM3D/branches/update_temporary/Makefile.in
===================================================================
--- seismo/3D/SPECFEM3D/branches/update_temporary/Makefile.in 2008-09-22 03:25:41 UTC (rev 12934)
+++ seismo/3D/SPECFEM3D/branches/update_temporary/Makefile.in 2008-09-22 15:01:36 UTC (rev 12935)
@@ -122,6 +122,7 @@
$O/save_arrays_solver.o \
$O/save_header_file.o \
$O/socal_model.o \
+ $O/sort_array_coordinates.o \
$O/utm_geo.o \
$O/write_AVS_DX_global_data.o \
$O/write_AVS_DX_global_faces_data.o \
@@ -420,6 +421,9 @@
$O/socal_model.o: constants.h socal_model.f90
${FCCOMPILE_CHECK} -c -o $O/socal_model.o socal_model.f90
+$O/sort_array_coordinates.o: constants.h sort_array_coordinates.f90
+ ${FCCOMPILE_CHECK} -c -o $O/sort_array_coordinates.o sort_array_coordinates.f90
+
$O/aniso_model.o: constants.h aniso_model.f90
${FCCOMPILE_CHECK} -c -o $O/aniso_model.o aniso_model.f90
Modified: seismo/3D/SPECFEM3D/branches/update_temporary/constants.h.in
===================================================================
--- seismo/3D/SPECFEM3D/branches/update_temporary/constants.h.in 2008-09-22 03:25:41 UTC (rev 12934)
+++ seismo/3D/SPECFEM3D/branches/update_temporary/constants.h.in 2008-09-22 15:01:36 UTC (rev 12935)
@@ -124,9 +124,18 @@
! deltat
integer, parameter :: DT_ext_mesh = 0.001d0
+
! NSTEP
double precision, parameter :: NSTEP_ext_mesh = 100
+! number of nodes per element as provided by the external mesh
+ integer, parameter :: ESIZE = 8
+
+! geometry tolerance parameter to calculate number of independent grid points
+! sensitive to actual size of model, assumes reference sphere of radius 1
+! this is an absolute value for normalized coordinates in the Earth
+ double precision, parameter :: SMALLVAL_TOL = 1.d-10
+
!------------------------------------------------------
!----------- do not modify anything below -------------
!------------------------------------------------------
Modified: seismo/3D/SPECFEM3D/branches/update_temporary/create_regions_mesh.f90
===================================================================
--- seismo/3D/SPECFEM3D/branches/update_temporary/create_regions_mesh.f90 2008-09-22 03:25:41 UTC (rev 12934)
+++ seismo/3D/SPECFEM3D/branches/update_temporary/create_regions_mesh.f90 2008-09-22 15:01:36 UTC (rev 12935)
@@ -1153,3 +1153,972 @@
end subroutine create_regions_mesh
+!
+!----
+!
+
+ subroutine create_regions_mesh_ext_mesh(ibool, &
+ xstore,ystore,zstore,npx,npy,iproc_xi,iproc_eta,nspec, &
+ volume_local,area_local_bottom,area_local_top, &
+ NGLOB_AB,npointot, &
+ myrank,LOCAL_PATH, &
+ nnodes_ext_mesh,nelmnts_ext_mesh, &
+ nodes_coords_ext_mesh,elmnts_ext_mesh,mat_ext_mesh, &
+ ninterface_ext_mesh,max_interface_size_ext_mesh, &
+ my_neighbours_ext_mesh,my_nelmnts_neighbours_ext_mesh,my_interfaces_ext_mesh, &
+ ibool_interfaces_ext_mesh,nibool_interfaces_ext_mesh &
+ )
+
+! create the different regions of the mesh
+
+ implicit none
+
+ include "constants.h"
+
+! number of spectral elements in each block
+ integer nspec
+
+ integer npx,npy
+ integer npointot
+
+ character(len=150) LOCAL_PATH
+
+! arrays with the mesh
+ double precision, dimension(NGLLX,NGLLY,NGLLZ,nspec) :: xstore,ystore,zstore
+
+ double precision xstore_local(NGLLX,NGLLY,NGLLZ)
+ double precision ystore_local(NGLLX,NGLLY,NGLLZ)
+ double precision zstore_local(NGLLX,NGLLY,NGLLZ)
+
+ double precision xmesh,ymesh,zmesh
+
+ integer ibool(NGLLX,NGLLY,NGLLZ,nspec)
+
+! data from the external mesh
+ integer :: nnodes_ext_mesh,nelmnts_ext_mesh
+ double precision, dimension(NDIM,nnodes_ext_mesh) :: nodes_coords_ext_mesh
+ integer, dimension(ESIZE,nelmnts_ext_mesh) :: elmnts_ext_mesh
+ integer, dimension(nelmnts_ext_mesh) :: mat_ext_mesh
+ double precision, external :: materials_ext_mesh
+ integer :: ninterface_ext_mesh,max_interface_size_ext_mesh
+ integer, dimension(ninterface_ext_mesh) :: my_neighbours_ext_mesh
+ integer, dimension(ninterface_ext_mesh) :: my_nelmnts_neighbours_ext_mesh
+ integer, dimension(6,max_interface_size_ext_mesh,ninterface_ext_mesh) :: my_interfaces_ext_mesh
+ integer, dimension(NGLLX*NGLLX*max_interface_size_ext_mesh,ninterface_ext_mesh) :: ibool_interfaces_ext_mesh
+ integer, dimension(ninterface_ext_mesh) :: nibool_interfaces_ext_mesh
+
+! for MPI buffers
+ integer, dimension(:), allocatable :: reorder_interface_ext_mesh,ind_ext_mesh,ninseg_ext_mesh,iwork_ext_mesh
+ integer, dimension(:), allocatable :: nibool_interfaces_ext_mesh_true
+ integer, dimension(:,:), allocatable :: ibool_interfaces_ext_mesh_dummy
+ integer, dimension(:), allocatable :: ibool_interface_ext_mesh_dummy
+ double precision, dimension(:), allocatable :: work_ext_mesh
+ real(kind=CUSTOM_REAL), dimension(:), allocatable :: xstore_dummy
+ real(kind=CUSTOM_REAL), dimension(:), allocatable :: ystore_dummy
+ real(kind=CUSTOM_REAL), dimension(:), allocatable :: zstore_dummy
+
+! Gauss-Lobatto-Legendre points and weights of integration
+ double precision, dimension(:), allocatable :: xigll,yigll,zigll,wxgll,wygll,wzgll
+
+! 3D shape functions and their derivatives
+ double precision, dimension(:,:,:,:), allocatable :: shape3D
+ double precision, dimension(:,:,:,:,:), allocatable :: dershape3D
+
+ double precision xelm(NGNOD)
+ double precision yelm(NGNOD)
+ double precision zelm(NGNOD)
+
+! the jacobian
+ real(kind=CUSTOM_REAL) jacobianl
+
+! arrays with mesh parameters
+ real(kind=CUSTOM_REAL), dimension(:,:,:,:), allocatable :: xixstore,xiystore,xizstore, &
+ etaxstore,etaystore,etazstore,gammaxstore,gammaystore,gammazstore,jacobianstore
+
+! for model density
+ real(kind=CUSTOM_REAL), dimension(:,:,:,:), allocatable :: rhostore,kappastore,mustore
+
+! proc numbers for MPI
+ integer myrank
+
+! check area and volume of the final mesh
+ double precision weight
+ double precision area_local_bottom,area_local_top
+ double precision volume_local
+
+! variables for creating array ibool (some arrays also used for AVS or DX files)
+ integer, dimension(:), allocatable :: iglob,locval
+ logical, dimension(:), allocatable :: ifseg
+ double precision, dimension(:), allocatable :: xp,yp,zp
+
+ integer nglob,NGLOB_AB
+ integer ieoff,ilocnum
+ integer ier
+ integer iinterface
+
+! mass matrix
+ real(kind=CUSTOM_REAL), dimension(:), allocatable :: rmass
+
+! ---------------------------
+
+! name of the database file
+ character(len=150) prname
+
+ integer i,j,k,ia,ispec,iglobnum,itype_element
+ integer iproc_xi,iproc_eta
+
+ double precision rho,vp,vs
+ double precision c11,c12,c13,c14,c15,c16,c22,c23,c24,c25,c26,c33,c34,c35,c36,c44,c45,c46,c55,c56,c66
+
+! for the harvard 3D salton sea model
+ double precision :: umesh, vmesh, wmesh, vp_st, vs_st, rho_st
+
+! mask to sort ibool
+ integer, dimension(:), allocatable :: mask_ibool
+ integer, dimension(:,:,:,:), allocatable :: copy_ibool_ori
+ integer :: inumber
+
+! **************
+
+! create the name for the database of the current slide and region
+ call create_name_database(prname,myrank,LOCAL_PATH)
+
+! Gauss-Lobatto-Legendre points of integration
+ allocate(xigll(NGLLX))
+ allocate(yigll(NGLLY))
+ allocate(zigll(NGLLZ))
+
+! Gauss-Lobatto-Legendre weights of integration
+ allocate(wxgll(NGLLX))
+ allocate(wygll(NGLLY))
+ allocate(wzgll(NGLLZ))
+
+! 3D shape functions and their derivatives
+ allocate(shape3D(NGNOD,NGLLX,NGLLY,NGLLZ))
+ allocate(dershape3D(NDIM,NGNOD,NGLLX,NGLLY,NGLLZ))
+
+! array with model density
+ allocate(rhostore(NGLLX,NGLLY,NGLLZ,nspec))
+ allocate(kappastore(NGLLX,NGLLY,NGLLZ,nspec))
+ allocate(mustore(NGLLX,NGLLY,NGLLZ,nspec))
+
+! arrays with mesh parameters
+ allocate(xixstore(NGLLX,NGLLY,NGLLZ,nspec))
+ allocate(xiystore(NGLLX,NGLLY,NGLLZ,nspec))
+ allocate(xizstore(NGLLX,NGLLY,NGLLZ,nspec))
+ allocate(etaxstore(NGLLX,NGLLY,NGLLZ,nspec))
+ allocate(etaystore(NGLLX,NGLLY,NGLLZ,nspec))
+ allocate(etazstore(NGLLX,NGLLY,NGLLZ,nspec))
+ allocate(gammaxstore(NGLLX,NGLLY,NGLLZ,nspec))
+ allocate(gammaystore(NGLLX,NGLLY,NGLLZ,nspec))
+ allocate(gammazstore(NGLLX,NGLLY,NGLLZ,nspec))
+ allocate(jacobianstore(NGLLX,NGLLY,NGLLZ,nspec))
+
+! set up coordinates of the Gauss-Lobatto-Legendre points
+ call zwgljd(xigll,wxgll,NGLLX,GAUSSALPHA,GAUSSBETA)
+ call zwgljd(yigll,wygll,NGLLY,GAUSSALPHA,GAUSSBETA)
+ call zwgljd(zigll,wzgll,NGLLZ,GAUSSALPHA,GAUSSBETA)
+
+! if number of points is odd, the middle abscissa is exactly zero
+ if(mod(NGLLX,2) /= 0) xigll((NGLLX-1)/2+1) = ZERO
+ if(mod(NGLLY,2) /= 0) yigll((NGLLY-1)/2+1) = ZERO
+ if(mod(NGLLZ,2) /= 0) zigll((NGLLZ-1)/2+1) = ZERO
+
+! get the 3-D shape functions
+ call get_shape3D(myrank,shape3D,dershape3D,xigll,yigll,zigll)
+
+! allocate memory for arrays
+ allocate(iglob(npointot))
+ allocate(locval(npointot))
+ allocate(ifseg(npointot))
+ allocate(xp(npointot))
+ allocate(yp(npointot))
+ allocate(zp(npointot))
+
+!---
+
+ xstore(:,:,:,:) = 0.d0
+ ystore(:,:,:,:) = 0.d0
+ zstore(:,:,:,:) = 0.d0
+
+ do ispec = 1, nspec
+ !call get_xyzelm(xelm, yelm, zelm, ispec, elmnts_ext_mesh, nodes_coords_ext_mesh, nspec, nnodes_ext_mesh)
+ do ia = 1,NGNOD
+ xelm(ia) = nodes_coords_ext_mesh(1,elmnts_ext_mesh(ia,ispec))
+ yelm(ia) = nodes_coords_ext_mesh(2,elmnts_ext_mesh(ia,ispec))
+ zelm(ia) = nodes_coords_ext_mesh(3,elmnts_ext_mesh(ia,ispec))
+ enddo
+
+ call calc_jacobian(myrank,xixstore,xiystore,xizstore, &
+ etaxstore,etaystore,etazstore, &
+ gammaxstore,gammaystore,gammazstore,jacobianstore, &
+ xstore,ystore,zstore, &
+ xelm,yelm,zelm,shape3D,dershape3D,ispec,nspec)
+
+ enddo
+
+! kappastore and mustore
+ do ispec = 1, nspec
+ do k = 1, NGLLZ
+ do j = 1, NGLLY
+ do i = 1, NGLLX
+ kappastore(i,j,k,ispec) = materials_ext_mesh(1,mat_ext_mesh(ispec))* &
+ (materials_ext_mesh(2,mat_ext_mesh(ispec))*materials_ext_mesh(2,mat_ext_mesh(ispec)) - &
+ 4.d0*materials_ext_mesh(3,mat_ext_mesh(ispec))*materials_ext_mesh(3,mat_ext_mesh(ispec))/3.d0)
+ mustore(i,j,k,ispec) = materials_ext_mesh(1,mat_ext_mesh(ispec))*materials_ext_mesh(3,mat_ext_mesh(ispec))*&
+ materials_ext_mesh(3,mat_ext_mesh(ispec))
+ enddo
+ enddo
+ enddo
+ enddo
+
+ locval = 0
+ ifseg = .false.
+ xp = 0.d0
+ yp = 0.d0
+ zp = 0.d0
+
+ do ispec=1,nspec
+ ieoff = NGLLX * NGLLY * NGLLZ * (ispec-1)
+ ilocnum = 0
+ do k=1,NGLLZ
+ do j=1,NGLLY
+ do i=1,NGLLX
+ ilocnum = ilocnum + 1
+ xp(ilocnum+ieoff) = xstore(i,j,k,ispec)
+ yp(ilocnum+ieoff) = ystore(i,j,k,ispec)
+ zp(ilocnum+ieoff) = zstore(i,j,k,ispec)
+ enddo
+ enddo
+ enddo
+ enddo
+
+ call get_global(nspec,xp,yp,zp,ibool,locval,ifseg,nglob,npointot, &
+ minval(nodes_coords_ext_mesh(1,:)),maxval(nodes_coords_ext_mesh(1,:)))
+
+ deallocate(xp,stat=ier); if(ier /= 0) stop 'error in deallocate'
+ deallocate(yp,stat=ier); if(ier /= 0) stop 'error in deallocate'
+ deallocate(zp,stat=ier); if(ier /= 0) stop 'error in deallocate'
+ deallocate(locval,stat=ier); if(ier /= 0) stop 'error in deallocate'
+ deallocate(ifseg,stat=ier); if(ier /= 0) stop 'error in deallocate'
+
+!
+!- we can create a new indirect addressing to reduce cache misses
+!
+ allocate(copy_ibool_ori(NGLLX,NGLLY,NGLLZ,nspec),stat=ier); if(ier /= 0) stop 'error in allocate'
+ allocate(mask_ibool(nglob),stat=ier); if(ier /= 0) stop 'error in allocate'
+
+ mask_ibool(:) = -1
+ copy_ibool_ori(:,:,:,:) = ibool(:,:,:,:)
+
+ inumber = 0
+ do ispec=1,nspec
+ do k=1,NGLLZ
+ do j=1,NGLLY
+ do i=1,NGLLX
+ if(mask_ibool(copy_ibool_ori(i,j,k,ispec)) == -1) then
+! create a new point
+ inumber = inumber + 1
+ ibool(i,j,k,ispec) = inumber
+ mask_ibool(copy_ibool_ori(i,j,k,ispec)) = inumber
+ else
+! use an existing point created previously
+ ibool(i,j,k,ispec) = mask_ibool(copy_ibool_ori(i,j,k,ispec))
+ endif
+ enddo
+ enddo
+ enddo
+ enddo
+
+ deallocate(copy_ibool_ori,stat=ier); if(ier /= 0) stop 'error in deallocate'
+ deallocate(mask_ibool,stat=ier); if(ier /= 0) stop 'error in deallocate'
+
+ allocate(xstore_dummy(nglob))
+ allocate(ystore_dummy(nglob))
+ allocate(zstore_dummy(nglob))
+ do ispec = 1, nspec
+ do k = 1, NGLLZ
+ do j = 1, NGLLY
+ do i = 1, NGLLX
+ iglobnum = ibool(i,j,k,ispec)
+ xstore_dummy(iglobnum) = xstore(i,j,k,ispec)
+ enddo
+ enddo
+ enddo
+ enddo
+
+ do ispec = 1, nspec
+ do k = 1, NGLLZ
+ do j = 1, NGLLY
+ do i = 1, NGLLX
+ iglobnum = ibool(i,j,k,ispec)
+ ystore_dummy(iglobnum) = ystore(i,j,k,ispec)
+ enddo
+ enddo
+ enddo
+ enddo
+
+ do ispec = 1, nspec
+ do k = 1, NGLLZ
+ do j = 1, NGLLY
+ do i = 1, NGLLX
+ iglobnum = ibool(i,j,k,ispec)
+ zstore_dummy(iglobnum) = zstore(i,j,k,ispec)
+ enddo
+ enddo
+ enddo
+ enddo
+
+! creating mass matrix (will be fully assembled with MPI in the solver)
+ allocate(rmass(nglob))
+ rmass(:) = 0._CUSTOM_REAL
+
+ do ispec=1,nspec
+ do k=1,NGLLZ
+ do j=1,NGLLY
+ do i=1,NGLLX
+ weight=wxgll(i)*wygll(j)*wzgll(k)
+ iglobnum=ibool(i,j,k,ispec)
+
+ jacobianl=jacobianstore(i,j,k,ispec)
+
+! distinguish between single and double precision for reals
+ if(CUSTOM_REAL == SIZE_REAL) then
+ rmass(iglobnum) = rmass(iglobnum) + &
+ sngl(dble(materials_ext_mesh(1,mat_ext_mesh(ispec))) * dble(jacobianl) * weight)
+ else
+ rmass(iglobnum) = rmass(iglobnum) + materials_ext_mesh(1,mat_ext_mesh(ispec)) * jacobianl * weight
+ endif
+
+ enddo
+ enddo
+ enddo
+ enddo
+
+ call prepare_assemble_MPI (nelmnts_ext_mesh,ibool, &
+ elmnts_ext_mesh, ESIZE, &
+ nglob, &
+ ninterface_ext_mesh, max_interface_size_ext_mesh, &
+ my_nelmnts_neighbours_ext_mesh, my_interfaces_ext_mesh, &
+ ibool_interfaces_ext_mesh, &
+ nibool_interfaces_ext_mesh &
+ )
+
+! sort ibool comm buffers lexicographically
+ allocate(nibool_interfaces_ext_mesh_true(ninterface_ext_mesh))
+
+ do iinterface = 1, ninterface_ext_mesh
+
+ allocate(xp(nibool_interfaces_ext_mesh(iinterface)))
+ allocate(yp(nibool_interfaces_ext_mesh(iinterface)))
+ allocate(zp(nibool_interfaces_ext_mesh(iinterface)))
+ allocate(locval(nibool_interfaces_ext_mesh(iinterface)))
+ allocate(ifseg(nibool_interfaces_ext_mesh(iinterface)))
+ allocate(reorder_interface_ext_mesh(nibool_interfaces_ext_mesh(iinterface)))
+ allocate(ibool_interface_ext_mesh_dummy(nibool_interfaces_ext_mesh(iinterface)))
+ allocate(ind_ext_mesh(nibool_interfaces_ext_mesh(iinterface)))
+ allocate(ninseg_ext_mesh(nibool_interfaces_ext_mesh(iinterface)))
+ allocate(iwork_ext_mesh(nibool_interfaces_ext_mesh(iinterface)))
+ allocate(work_ext_mesh(nibool_interfaces_ext_mesh(iinterface)))
+
+ do ilocnum = 1, nibool_interfaces_ext_mesh(iinterface)
+ xp(ilocnum) = xstore_dummy(ibool_interfaces_ext_mesh(ilocnum,iinterface))
+ yp(ilocnum) = ystore_dummy(ibool_interfaces_ext_mesh(ilocnum,iinterface))
+ zp(ilocnum) = zstore_dummy(ibool_interfaces_ext_mesh(ilocnum,iinterface))
+ enddo
+
+ call sort_array_coordinates(nibool_interfaces_ext_mesh(iinterface),xp,yp,zp, &
+ ibool_interfaces_ext_mesh(1:nibool_interfaces_ext_mesh(iinterface),iinterface), &
+ reorder_interface_ext_mesh,locval,ifseg,nibool_interfaces_ext_mesh_true(iinterface), &
+ ind_ext_mesh,ninseg_ext_mesh,iwork_ext_mesh,work_ext_mesh)
+
+ deallocate(xp)
+ deallocate(yp)
+ deallocate(zp)
+ deallocate(locval)
+ deallocate(ifseg)
+ deallocate(reorder_interface_ext_mesh)
+ deallocate(ibool_interface_ext_mesh_dummy)
+ deallocate(ind_ext_mesh)
+ deallocate(ninseg_ext_mesh)
+ deallocate(iwork_ext_mesh)
+ deallocate(work_ext_mesh)
+
+ enddo
+
+! save the binary files
+ call create_name_database(prname,myrank,LOCAL_PATH)
+ open(unit=IOUT,file=prname(1:len_trim(prname))//'external_mesh.bin',status='unknown',action='write',form='unformatted')
+ write(IOUT) nspec
+ write(IOUT) nglob
+
+ write(IOUT) xixstore
+ write(IOUT) xiystore
+ write(IOUT) xizstore
+ write(IOUT) etaxstore
+ write(IOUT) etaystore
+ write(IOUT) etazstore
+ write(IOUT) gammaxstore
+ write(IOUT) gammaystore
+ write(IOUT) gammazstore
+
+ write(IOUT) jacobianstore
+ write(IOUT) kappastore
+ write(IOUT) mustore
+
+ write(IOUT) rmass
+
+ write(IOUT) ibool
+
+ write(IOUT) xstore_dummy
+ write(IOUT) ystore_dummy
+ write(IOUT) zstore_dummy
+
+ write(IOUT) ninterface_ext_mesh
+ write(IOUT) maxval(nibool_interfaces_ext_mesh)
+ write(IOUT) my_neighbours_ext_mesh
+ write(IOUT) nibool_interfaces_ext_mesh
+ allocate(ibool_interfaces_ext_mesh_dummy(maxval(nibool_interfaces_ext_mesh),ninterface_ext_mesh))
+ do i = 1, ninterface_ext_mesh
+ ibool_interfaces_ext_mesh_dummy = ibool_interfaces_ext_mesh(1:maxval(nibool_interfaces_ext_mesh),:)
+ enddo
+ write(IOUT) ibool_interfaces_ext_mesh_dummy
+ close(IOUT)
+
+ end subroutine create_regions_mesh_ext_mesh
+
+!
+!----
+!
+
+ double precision function materials_ext_mesh(i,j)
+
+ implicit none
+
+ integer :: i,j
+
+ select case (j)
+ case (1)
+ select case (i)
+ case (1)
+ materials_ext_mesh = 2000.d0
+ case (2)
+ materials_ext_mesh = 3000.d0
+ case (3)
+ materials_ext_mesh = 1732.051d0
+ case default
+ call stop_all()
+ end select
+ case (2)
+ select case (i)
+ case (1)
+ materials_ext_mesh = 2000.d0
+ case (2)
+ materials_ext_mesh = 900.d0
+ case (3)
+ materials_ext_mesh = 500.d0
+ case default
+ call stop_all()
+ end select
+ case default
+ call stop_all()
+ end select
+
+ end function materials_ext_mesh
+
+!
+!----
+!
+
+subroutine prepare_assemble_MPI (nelmnts,ibool, &
+ knods, ngnode, &
+ npoin, &
+ ninterface, max_interface_size, &
+ my_nelmnts_neighbours, my_interfaces, &
+ ibool_interfaces_asteroid, &
+ nibool_interfaces_asteroid &
+ )
+
+ implicit none
+
+ include 'constants.h'
+
+ integer, intent(in) :: nelmnts, npoin, ngnode
+ integer, dimension(ngnode,nelmnts), intent(in) :: knods
+ integer, dimension(NGLLX,NGLLY,NGLLZ,nelmnts), intent(in) :: ibool
+
+ integer :: ninterface
+ integer :: max_interface_size
+ integer, dimension(ninterface) :: my_nelmnts_neighbours
+ integer, dimension(6,max_interface_size,ninterface) :: my_interfaces
+ integer, dimension(NGLLX*NGLLX*max_interface_size,ninterface) :: &
+ ibool_interfaces_asteroid
+ integer, dimension(ninterface) :: &
+ nibool_interfaces_asteroid
+
+ integer :: num_interface
+ integer :: ispec_interface
+
+ logical, dimension(npoin) :: mask_ibool_asteroid
+
+ integer :: ixmin, ixmax
+ integer :: iymin, iymax
+ integer :: izmin, izmax
+ integer, dimension(ngnode) :: n
+ integer :: e1, e2, e3, e4
+ integer :: type
+ integer :: ispec
+
+ integer :: k
+ integer :: npoin_interface_asteroid
+
+ integer :: ix,iy,iz
+
+
+ ibool_interfaces_asteroid(:,:) = 0
+ nibool_interfaces_asteroid(:) = 0
+
+ do num_interface = 1, ninterface
+ npoin_interface_asteroid = 0
+ mask_ibool_asteroid(:) = .false.
+
+ do ispec_interface = 1, my_nelmnts_neighbours(num_interface)
+ ispec = my_interfaces(1,ispec_interface,num_interface)
+ type = my_interfaces(2,ispec_interface,num_interface)
+ do k = 1, ngnode
+ n(k) = knods(k,ispec)
+ end do
+ e1 = my_interfaces(3,ispec_interface,num_interface)
+ e2 = my_interfaces(4,ispec_interface,num_interface)
+ e3 = my_interfaces(5,ispec_interface,num_interface)
+ e4 = my_interfaces(6,ispec_interface,num_interface)
+ call get_edge(ngnode, n, type, e1, e2, e3, e4, ixmin, ixmax, iymin, iymax, izmin, izmax)
+
+ do iz = min(izmin,izmax), max(izmin,izmax)
+ do iy = min(iymin,iymax), max(iymin,iymax)
+ do ix = min(ixmin,ixmax), max(ixmin,ixmax)
+
+ if(.not. mask_ibool_asteroid(ibool(ix,iy,iz,ispec))) then
+ mask_ibool_asteroid(ibool(ix,iy,iz,ispec)) = .true.
+ npoin_interface_asteroid = npoin_interface_asteroid + 1
+ ibool_interfaces_asteroid(npoin_interface_asteroid,num_interface)=&
+ ibool(ix,iy,iz,ispec)
+ end if
+ end do
+ end do
+ end do
+
+ end do
+ nibool_interfaces_asteroid(num_interface) = npoin_interface_asteroid
+
+
+ end do
+
+end subroutine prepare_assemble_MPI
+
+!
+!----
+!
+
+subroutine get_edge ( ngnode, n, type, e1, e2, e3, e4, ixmin, ixmax, iymin, iymax, izmin, izmax )
+
+ implicit none
+
+ include "constants.h"
+
+ integer, intent(in) :: ngnode
+ integer, dimension(ngnode), intent(in) :: n
+ integer, intent(in) :: type, e1, e2, e3, e4
+ integer, intent(out) :: ixmin, ixmax, iymin, iymax, izmin, izmax
+
+ integer, dimension(4) :: en
+ integer :: valence, i
+
+ if ( type == 1 ) then
+ if ( e1 == n(1) ) then
+ ixmin = 1
+ ixmax = 1
+ iymin = 1
+ iymax = 1
+ izmin = 1
+ izmax = 1
+ end if
+ if ( e1 == n(2) ) then
+ ixmin = NGLLX
+ ixmax = NGLLX
+ iymin = 1
+ iymax = 1
+ izmin = 1
+ izmax = 1
+ end if
+ if ( e1 == n(3) ) then
+ ixmin = NGLLX
+ ixmax = NGLLX
+ iymin = NGLLY
+ iymax = NGLLY
+ izmin = 1
+ izmax = 1
+ end if
+ if ( e1 == n(4) ) then
+ ixmin = 1
+ ixmax = 1
+ iymin = NGLLY
+ iymax = NGLLY
+ izmin = 1
+ izmax = 1
+ end if
+ if ( e1 == n(5) ) then
+ ixmin = 1
+ ixmax = 1
+ iymin = 1
+ iymax = 1
+ izmin = NGLLZ
+ izmax = NGLLZ
+ end if
+ if ( e1 == n(6) ) then
+ ixmin = NGLLX
+ ixmax = NGLLX
+ iymin = 1
+ iymax = 1
+ izmin = NGLLZ
+ izmax = NGLLZ
+ end if
+ if ( e1 == n(7) ) then
+ ixmin = NGLLX
+ ixmax = NGLLX
+ iymin = NGLLY
+ iymax = NGLLY
+ izmin = NGLLZ
+ izmax = NGLLZ
+ end if
+ if ( e1 == n(8) ) then
+ ixmin = 1
+ ixmax = 1
+ iymin = NGLLY
+ iymax = NGLLY
+ izmin = NGLLZ
+ izmax = NGLLZ
+ end if
+ else
+ if ( type == 2 ) then
+ if ( e1 == n(1) ) then
+ ixmin = 1
+ iymin = 1
+ izmin = 1
+ if ( e2 == n(2) ) then
+ ixmax = NGLLX
+ iymax = 1
+ izmax = 1
+ end if
+ if ( e2 == n(4) ) then
+ ixmax = 1
+ iymax = NGLLY
+ izmax = 1
+ end if
+ if ( e2 == n(5) ) then
+ ixmax = 1
+ iymax = 1
+ izmax = NGLLZ
+ end if
+ end if
+ if ( e1 == n(2) ) then
+ ixmin = NGLLX
+ iymin = 1
+ izmin = 1
+ if ( e2 == n(3) ) then
+ ixmax = NGLLX
+ iymax = NGLLY
+ izmax = 1
+ end if
+ if ( e2 == n(1) ) then
+ ixmax = 1
+ iymax = 1
+ izmax = 1
+ end if
+ if ( e2 == n(6) ) then
+ ixmax = NGLLX
+ iymax = 1
+ izmax = NGLLZ
+ end if
+
+ end if
+ if ( e1 == n(3) ) then
+ ixmin = NGLLX
+ iymin = NGLLY
+ izmin = 1
+ if ( e2 == n(4) ) then
+ ixmax = 1
+ iymax = NGLLY
+ izmax = 1
+ end if
+ if ( e2 == n(2) ) then
+ ixmax = NGLLX
+ iymax = 1
+ izmax = 1
+ end if
+ if ( e2 == n(7) ) then
+ ixmax = NGLLX
+ iymax = NGLLY
+ izmax = NGLLZ
+ end if
+ end if
+ if ( e1 == n(4) ) then
+ ixmin = 1
+ iymin = NGLLY
+ izmin = 1
+ if ( e2 == n(1) ) then
+ ixmax = 1
+ iymax = 1
+ izmax = 1
+ end if
+ if ( e2 == n(3) ) then
+ ixmax = NGLLX
+ iymax = NGLLY
+ izmax = 1
+ end if
+ if ( e2 == n(8) ) then
+ ixmax = 1
+ iymax = NGLLY
+ izmax = NGLLZ
+ end if
+ end if
+ if ( e1 == n(5) ) then
+ ixmin = 1
+ iymin = 1
+ izmin = NGLLZ
+ if ( e2 == n(1) ) then
+ ixmax = 1
+ iymax = 1
+ izmax = 1
+ end if
+ if ( e2 == n(6) ) then
+ ixmax = NGLLX
+ iymax = 1
+ izmax = NGLLZ
+ end if
+ if ( e2 == n(8) ) then
+ ixmax = 1
+ iymax = NGLLY
+ izmax = NGLLZ
+ end if
+ end if
+ if ( e1 == n(6) ) then
+ ixmin = NGLLX
+ iymin = 1
+ izmin = NGLLZ
+ if ( e2 == n(2) ) then
+ ixmax = NGLLX
+ iymax = 1
+ izmax = 1
+ end if
+ if ( e2 == n(7) ) then
+ ixmax = NGLLX
+ iymax = NGLLY
+ izmax = NGLLZ
+ end if
+ if ( e2 == n(5) ) then
+ ixmax = 1
+ iymax = 1
+ izmax = NGLLZ
+ end if
+ end if
+ if ( e1 == n(7) ) then
+ ixmin = NGLLX
+ iymin = NGLLY
+ izmin = NGLLZ
+ if ( e2 == n(3) ) then
+ ixmax = NGLLX
+ iymax = NGLLY
+ izmax = 1
+ end if
+ if ( e2 == n(8) ) then
+ ixmax = 1
+ iymax = NGLLY
+ izmax = NGLLZ
+ end if
+ if ( e2 == n(6) ) then
+ ixmax = NGLLX
+ iymax = 1
+ izmax = NGLLZ
+ end if
+ end if
+ if ( e1 == n(8) ) then
+ ixmin = 1
+ iymin = NGLLY
+ izmin = NGLLZ
+ if ( e2 == n(4) ) then
+ ixmax = 1
+ iymax = NGLLY
+ izmax = 1
+ end if
+ if ( e2 == n(5) ) then
+ ixmax = 1
+ iymax = 1
+ izmax = NGLLZ
+ end if
+ if ( e2 == n(7) ) then
+ ixmax = NGLLX
+ iymax = NGLLY
+ izmax = NGLLZ
+ end if
+ end if
+
+ else
+ if (type == 4) then
+ en(1) = e1
+ en(2) = e2
+ en(3) = e3
+ en(4) = e4
+
+ valence = 0
+ do i = 1, 4
+ if ( en(i) == n(1)) then
+ valence = valence+1
+ endif
+ if ( en(i) == n(2)) then
+ valence = valence+1
+ endif
+ if ( en(i) == n(3)) then
+ valence = valence+1
+ endif
+ if ( en(i) == n(4)) then
+ valence = valence+1
+ endif
+ enddo
+ if ( valence == 4 ) then
+ ixmin = 1
+ iymin = 1
+ izmin = 1
+ ixmax = NGLLX
+ iymax = NGLLY
+ izmax = 1
+ endif
+
+ valence = 0
+ do i = 1, 4
+ if ( en(i) == n(1)) then
+ valence = valence+1
+ endif
+ if ( en(i) == n(2)) then
+ valence = valence+1
+ endif
+ if ( en(i) == n(5)) then
+ valence = valence+1
+ endif
+ if ( en(i) == n(6)) then
+ valence = valence+1
+ endif
+ enddo
+ if ( valence == 4 ) then
+ ixmin = 1
+ iymin = 1
+ izmin = 1
+ ixmax = NGLLX
+ iymax = 1
+ izmax = NGLLZ
+ endif
+
+ valence = 0
+ do i = 1, 4
+ if ( en(i) == n(2)) then
+ valence = valence+1
+ endif
+ if ( en(i) == n(3)) then
+ valence = valence+1
+ endif
+ if ( en(i) == n(6)) then
+ valence = valence+1
+ endif
+ if ( en(i) == n(7)) then
+ valence = valence+1
+ endif
+ enddo
+ if ( valence == 4 ) then
+ ixmin = NGLLX
+ iymin = 1
+ izmin = 1
+ ixmax = NGLLX
+ iymax = NGLLZ
+ izmax = NGLLZ
+ endif
+
+ valence = 0
+ do i = 1, 4
+ if ( en(i) == n(3)) then
+ valence = valence+1
+ endif
+ if ( en(i) == n(4)) then
+ valence = valence+1
+ endif
+ if ( en(i) == n(7)) then
+ valence = valence+1
+ endif
+ if ( en(i) == n(8)) then
+ valence = valence+1
+ endif
+ enddo
+ if ( valence == 4 ) then
+ ixmin = 1
+ iymin = NGLLY
+ izmin = 1
+ ixmax = NGLLX
+ iymax = NGLLY
+ izmax = NGLLZ
+ endif
+
+ valence = 0
+ do i = 1, 4
+ if ( en(i) == n(1)) then
+ valence = valence+1
+ endif
+ if ( en(i) == n(4)) then
+ valence = valence+1
+ endif
+ if ( en(i) == n(5)) then
+ valence = valence+1
+ endif
+ if ( en(i) == n(8)) then
+ valence = valence+1
+ endif
+ enddo
+ if ( valence == 4 ) then
+ ixmin = 1
+ iymin = 1
+ izmin = 1
+ ixmax = 1
+ iymax = NGLLY
+ izmax = NGLLZ
+ endif
+
+ valence = 0
+ do i = 1, 4
+ if ( en(i) == n(5)) then
+ valence = valence+1
+ endif
+ if ( en(i) == n(6)) then
+ valence = valence+1
+ endif
+ if ( en(i) == n(7)) then
+ valence = valence+1
+ endif
+ if ( en(i) == n(8)) then
+ valence = valence+1
+ endif
+ enddo
+ if ( valence == 4 ) then
+ ixmin = 1
+ iymin = 1
+ izmin = NGLLZ
+ ixmax = NGLLX
+ iymax = NGLLY
+ izmax = NGLLZ
+ endif
+
+ else
+ stop 'ERROR get_edge'
+ endif
+
+ end if
+ end if
+
+end subroutine get_edge
Modified: seismo/3D/SPECFEM3D/branches/update_temporary/meshfem3D.f90
===================================================================
--- seismo/3D/SPECFEM3D/branches/update_temporary/meshfem3D.f90 2008-09-22 03:25:41 UTC (rev 12934)
+++ seismo/3D/SPECFEM3D/branches/update_temporary/meshfem3D.f90 2008-09-22 15:01:36 UTC (rev 12935)
@@ -228,6 +228,25 @@
character(len=MAX_LENGTH_NETWORK_NAME) network_name
character(len=150) rec_filename,filtered_rec_filename,dummystring
+! for Databases of external meshes
+ character(len=150) prname
+ integer :: dummy_node
+ integer :: dummy_elmnt
+ integer :: ispec, inode, num_interface, ie
+ integer :: nnodes_ext_mesh, nelmnts_ext_mesh
+ integer :: ninterface_ext_mesh
+ integer :: max_interface_size_ext_mesh
+ integer, dimension(:), allocatable :: my_neighbours_ext_mesh
+ integer, dimension(:), allocatable :: my_nelmnts_neighbours_ext_mesh
+ integer, dimension(:,:,:), allocatable :: my_interfaces_ext_mesh
+ integer, dimension(:,:), allocatable :: ibool_interfaces_ext_mesh
+ integer, dimension(:,:), allocatable :: ibool_interfaces_ext_mesh_dummy
+ integer, dimension(:), allocatable :: ibool_interface_ext_mesh_dummy
+ integer, dimension(:), allocatable :: nibool_interfaces_ext_mesh
+ double precision, dimension(:,:), allocatable :: nodes_coords_ext_mesh
+ integer, dimension(:,:), allocatable :: elmnts_ext_mesh
+ integer, dimension(:), allocatable :: mat_ext_mesh
+
! ************** PROGRAM STARTS HERE **************
! sizeprocs returns number of processes started (should be equal to NPROC).
@@ -282,9 +301,16 @@
NSPEC2DMAX_XMIN_XMAX,NSPEC2DMAX_YMIN_YMAX,NSPEC2D_BOTTOM,NSPEC2D_TOP, &
NPOIN2DMAX_XMIN_XMAX,NPOIN2DMAX_YMIN_YMAX,NGLOB_AB,USE_REGULAR_MESH)
+! info about external mesh simulation
+! nlegoff -- should be put in compute_parameters and read_parameter_file for clarity
+ if (USE_EXTERNAL_MESH) then
+ NPROC = sizeprocs
+ endif
+
! check that the code is running with the requested nb of processes
if(sizeprocs /= NPROC) call exit_MPI(myrank,'wrong number of MPI processes')
+ if (.not. USE_EXTERNAL_MESH) then
! dynamic allocation of mesh arrays
allocate(rns(0:2*NER))
@@ -330,6 +356,8 @@
write(IMAIN,'(a1)',advance='yes') ' '
enddo
endif
+
+ endif ! end of (.not. USE_EXTERNAL_MESH)
if(myrank == 0) then
write(IMAIN,*) 'This is process ',myrank
@@ -362,6 +390,7 @@
! for the number of standard linear solids for attenuation
if(N_SLS /= 3) call exit_MPI(myrank,'number of SLS must be 3')
+ if (.not. USE_EXTERNAL_MESH) then
! check that Poisson's ratio in Gocad block is fine
if(VP_VS_RATIO_GOCAD_TOP < sqrt(2.) .or. VP_VS_RATIO_GOCAD_BOTTOM < sqrt(2.))&
call exit_MPI(myrank,'vp/vs ratio in Gocad block is too small')
@@ -383,6 +412,8 @@
if(mod(NEX_ETA/8,NPROC_ETA) /= 0) call exit_MPI(myrank,'NEX_ETA must be a multiple of 8*NPROC_ETA for a non-regular mesh')
endif
+ endif ! end of (.not. USE_EXTERNAL_MESH)
+
if(myrank == 0) then
write(IMAIN,*) 'region selected:'
@@ -498,6 +529,8 @@
close(55)
endif
+ if (.not. USE_EXTERNAL_MESH) then
+
! get addressing for this process
iproc_xi = iproc_xi_slice(myrank)
iproc_eta = iproc_eta_slice(myrank)
@@ -672,6 +705,8 @@
enddo
enddo
+ endif ! end of (.not. USE_EXTERNAL_MESH)
+
if(myrank == 0) then
write(IMAIN,*)
write(IMAIN,*) '**************************'
@@ -685,6 +720,46 @@
area_local_bottom = ZERO
area_local_top = ZERO
+! read databases about external mesh simulation
+! nlegoff --
+ if (USE_EXTERNAL_MESH) then
+ call create_name_database(prname,myrank,LOCAL_PATH)
+ open(unit=IIN,file=prname(1:len_trim(prname))//'Database',status='old',action='read',form='formatted')
+ read(IIN,*) nnodes_ext_mesh
+ allocate(nodes_coords_ext_mesh(NDIM,nnodes_ext_mesh))
+ do inode = 1, nnodes_ext_mesh
+ read(IIN,*) dummy_node, nodes_coords_ext_mesh(1,inode), nodes_coords_ext_mesh(2,inode), nodes_coords_ext_mesh(3,inode)
+ enddo
+
+ read(IIN,*) nelmnts_ext_mesh
+ allocate(elmnts_ext_mesh(esize,nelmnts_ext_mesh))
+ allocate(mat_ext_mesh(nelmnts_ext_mesh))
+ do ispec = 1, nelmnts_ext_mesh
+ read(IIN,*) dummy_elmnt, mat_ext_mesh(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)
+ enddo
+ NSPEC_AB = nelmnts_ext_mesh
+
+ read(IIN,*) ninterface_ext_mesh, max_interface_size_ext_mesh
+ allocate(my_neighbours_ext_mesh(ninterface_ext_mesh))
+ allocate(my_nelmnts_neighbours_ext_mesh(ninterface_ext_mesh))
+ allocate(my_interfaces_ext_mesh(6,max_interface_size_ext_mesh,ninterface_ext_mesh))
+ allocate(ibool_interfaces_ext_mesh(NGLLX*NGLLX*max_interface_size_ext_mesh,ninterface_ext_mesh))
+ allocate(nibool_interfaces_ext_mesh(ninterface_ext_mesh))
+ do num_interface = 1, ninterface_ext_mesh
+ read(IIN,*) my_neighbours_ext_mesh(num_interface), my_nelmnts_neighbours_ext_mesh(num_interface)
+ do ie = 1, my_nelmnts_neighbours_ext_mesh(num_interface)
+ read(IIN,*) my_interfaces_ext_mesh(1,ie,num_interface), my_interfaces_ext_mesh(2,ie,num_interface), &
+ my_interfaces_ext_mesh(3,ie,num_interface), my_interfaces_ext_mesh(4,ie,num_interface), &
+ my_interfaces_ext_mesh(5,ie,num_interface), my_interfaces_ext_mesh(6,ie,num_interface)
+ enddo
+ enddo
+
+ close(IIN)
+
+ endif
+
! assign theoretical number of elements
nspec = NSPEC_AB
@@ -705,6 +780,20 @@
if(ier /= 0) call exit_MPI(myrank,'not enough memory to allocate arrays')
! create all the regions of the mesh
+ if (USE_EXTERNAL_MESH) then
+ call create_regions_mesh_ext_mesh(ibool, &
+ xstore,ystore,zstore,npx,npy,iproc_xi,iproc_eta,nspec, &
+ volume_local,area_local_bottom,area_local_top, &
+ NGLOB_AB,npointot, &
+ myrank,LOCAL_PATH, &
+ nnodes_ext_mesh,nelmnts_ext_mesh, &
+ nodes_coords_ext_mesh,elmnts_ext_mesh,mat_ext_mesh, &
+ ninterface_ext_mesh,max_interface_size_ext_mesh, &
+ my_neighbours_ext_mesh,my_nelmnts_neighbours_ext_mesh,my_interfaces_ext_mesh, &
+ ibool_interfaces_ext_mesh,nibool_interfaces_ext_mesh &
+ )
+
+ else
call create_regions_mesh(xgrid,ygrid,zgrid,ibool,idoubling, &
xstore,ystore,zstore,npx,npy, &
iproc_xi,iproc_eta,nspec, &
@@ -722,7 +811,7 @@
IMPOSE_MINIMUM_VP_GOCAD,THICKNESS_TAPER_BLOCK_HR,THICKNESS_TAPER_BLOCK_MR,MOHO_MAP_LUPEI, &
ANISOTROPY,SAVE_MESH_FILES,SUPPRESS_UTM_PROJECTION, &
ORIG_LAT_TOPO,ORIG_LONG_TOPO,DEGREES_PER_CELL_TOPO,NX_TOPO,NY_TOPO,USE_REGULAR_MESH)
-
+ endif
! print min and max of topography included
if(TOPOGRAPHY) then
Added: seismo/3D/SPECFEM3D/branches/update_temporary/sort_array_coordinates.f90
===================================================================
--- seismo/3D/SPECFEM3D/branches/update_temporary/sort_array_coordinates.f90 (rev 0)
+++ seismo/3D/SPECFEM3D/branches/update_temporary/sort_array_coordinates.f90 2008-09-22 15:01:36 UTC (rev 12935)
@@ -0,0 +1,235 @@
+!=====================================================================
+!
+! S p e c f e m 3 D G l o b e V e r s i o n 4 . 0
+! --------------------------------------------------
+!
+! Main authors: Dimitri Komatitsch and Jeroen Tromp
+! Seismological Laboratory, California Institute of Technology, USA
+! and University of Pau / CNRS / INRIA, France
+! (c) California Institute of Technology and University of Pau / CNRS / INRIA
+! February 2008
+!
+! This program is free software; you can redistribute it and/or modify
+! it under the terms of the GNU General Public License as published by
+! the Free Software Foundation; either version 2 of the License, or
+! (at your option) any later version.
+!
+! This program is distributed in the hope that it will be useful,
+! but WITHOUT ANY WARRANTY; without even the implied warranty of
+! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+! GNU General Public License for more details.
+!
+! You should have received a copy of the GNU General Public License along
+! with this program; if not, write to the Free Software Foundation, Inc.,
+! 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
+!
+!=====================================================================
+
+! subroutines to sort MPI buffers to assemble between chunks
+
+ subroutine sort_array_coordinates(npointot,x,y,z,ibool,iglob,loc,ifseg,nglob,ind,ninseg,iwork,work)
+
+! this routine MUST be in double precision to avoid sensitivity
+! to roundoff errors in the coordinates of the points
+
+ implicit none
+
+ include "constants.h"
+
+ integer npointot,nglob
+
+ integer ibool(npointot),iglob(npointot),loc(npointot)
+ integer ind(npointot),ninseg(npointot)
+ logical ifseg(npointot)
+ double precision x(npointot),y(npointot),z(npointot)
+ integer iwork(npointot)
+ double precision work(npointot)
+
+ integer ipoin,i,j
+ integer nseg,ioff,iseg,ig
+ double precision xtol
+
+! establish initial pointers
+ do ipoin=1,npointot
+ loc(ipoin)=ipoin
+ enddo
+
+! define a tolerance, normalized radius is 1., so let's use a small value
+ xtol = SMALLVAL_TOL
+
+ ifseg(:)=.false.
+
+ nseg=1
+ ifseg(1)=.true.
+ ninseg(1)=npointot
+
+ do j=1,NDIM
+
+! sort within each segment
+ ioff=1
+ do iseg=1,nseg
+ if(j == 1) then
+
+ call rank_buffers(x(ioff),ind,ninseg(iseg))
+
+ else if(j == 2) then
+
+ call rank_buffers(y(ioff),ind,ninseg(iseg))
+
+ else
+
+ call rank_buffers(z(ioff),ind,ninseg(iseg))
+
+ endif
+
+ call swap_all_buffers(ibool(ioff),loc(ioff), &
+ x(ioff),y(ioff),z(ioff),iwork,work,ind,ninseg(iseg))
+
+ ioff=ioff+ninseg(iseg)
+ enddo
+
+! check for jumps in current coordinate
+ if(j == 1) then
+ do i=2,npointot
+ if(dabs(x(i)-x(i-1)) > xtol) ifseg(i)=.true.
+ enddo
+ else if(j == 2) then
+ do i=2,npointot
+ if(dabs(y(i)-y(i-1)) > xtol) ifseg(i)=.true.
+ enddo
+ else
+ do i=2,npointot
+ if(dabs(z(i)-z(i-1)) > xtol) ifseg(i)=.true.
+ enddo
+ endif
+
+! count up number of different segments
+ nseg=0
+ do i=1,npointot
+ if(ifseg(i)) then
+ nseg=nseg+1
+ ninseg(nseg)=1
+ else
+ ninseg(nseg)=ninseg(nseg)+1
+ endif
+ enddo
+ enddo
+
+! assign global node numbers (now sorted lexicographically)
+ ig=0
+ do i=1,npointot
+ if(ifseg(i)) ig=ig+1
+ iglob(loc(i))=ig
+ enddo
+
+ nglob=ig
+
+ end subroutine sort_array_coordinates
+
+! -------------------- library for sorting routine ------------------
+
+! sorting routines put here in same file to allow for inlining
+
+ subroutine rank_buffers(A,IND,N)
+!
+! Use Heap Sort (Numerical Recipes)
+!
+ implicit none
+
+ integer n
+ double precision A(n)
+ integer IND(n)
+
+ integer i,j,l,ir,indx
+ double precision q
+
+ do j=1,n
+ IND(j)=j
+ enddo
+
+ if(n == 1) return
+
+ L=n/2+1
+ ir=n
+ 100 CONTINUE
+ IF(l>1) THEN
+ l=l-1
+ indx=ind(l)
+ q=a(indx)
+ ELSE
+ indx=ind(ir)
+ q=a(indx)
+ ind(ir)=ind(1)
+ ir=ir-1
+ if (ir == 1) then
+ ind(1)=indx
+ return
+ endif
+ ENDIF
+ i=l
+ j=l+l
+ 200 CONTINUE
+ IF(J <= IR) THEN
+ IF(J < IR) THEN
+ IF(A(IND(j)) < A(IND(j+1))) j=j+1
+ ENDIF
+ IF (q < A(IND(j))) THEN
+ IND(I)=IND(J)
+ I=J
+ J=J+J
+ ELSE
+ J=IR+1
+ ENDIF
+ goto 200
+ ENDIF
+ IND(I)=INDX
+ goto 100
+ end subroutine rank_buffers
+
+! -------------------------------------------------------------------
+
+ subroutine swap_all_buffers(IA,IB,A,B,C,IW,W,ind,n)
+!
+! swap arrays IA, IB, A, B and C according to addressing in array IND
+!
+ implicit none
+
+ integer n
+
+ integer IND(n)
+ integer IA(n),IB(n),IW(n)
+ double precision A(n),B(n),C(n),W(n)
+
+ integer i
+
+ do i=1,n
+ W(i)=A(i)
+ IW(i)=IA(i)
+ enddo
+
+ do i=1,n
+ A(i)=W(ind(i))
+ IA(i)=IW(ind(i))
+ enddo
+
+ do i=1,n
+ W(i)=B(i)
+ IW(i)=IB(i)
+ enddo
+
+ do i=1,n
+ B(i)=W(ind(i))
+ IB(i)=IW(ind(i))
+ enddo
+
+ do i=1,n
+ W(i)=C(i)
+ enddo
+
+ do i=1,n
+ C(i)=W(ind(i))
+ enddo
+
+ end subroutine swap_all_buffers
+
+
More information about the cig-commits
mailing list