[cig-commits] r18311 - in seismo/3D/FAULT_SOURCE/branches: . decompose_mesh_SCOTCH
percygalvez at geodynamics.org
percygalvez at geodynamics.org
Tue May 3 18:12:20 PDT 2011
Author: percygalvez
Date: 2011-05-03 18:12:19 -0700 (Tue, 03 May 2011)
New Revision: 18311
Added:
seismo/3D/FAULT_SOURCE/branches/decompose_mesh_SCOTCH/
seismo/3D/FAULT_SOURCE/branches/decompose_mesh_SCOTCH/Makefile
seismo/3D/FAULT_SOURCE/branches/decompose_mesh_SCOTCH/README
seismo/3D/FAULT_SOURCE/branches/decompose_mesh_SCOTCH/compile_all.bash
seismo/3D/FAULT_SOURCE/branches/decompose_mesh_SCOTCH/decompose_mesh_SCOTCH.f90
seismo/3D/FAULT_SOURCE/branches/decompose_mesh_SCOTCH/fault_scotch.f90
seismo/3D/FAULT_SOURCE/branches/decompose_mesh_SCOTCH/part_decompose_mesh_SCOTCH.f90
seismo/3D/FAULT_SOURCE/branches/decompose_mesh_SCOTCH/program_decompose_mesh_SCOTCH.f90
seismo/3D/FAULT_SOURCE/branches/decompose_mesh_SCOTCH/scotchf.h
Log:
branches added
Added: seismo/3D/FAULT_SOURCE/branches/decompose_mesh_SCOTCH/Makefile
===================================================================
--- seismo/3D/FAULT_SOURCE/branches/decompose_mesh_SCOTCH/Makefile (rev 0)
+++ seismo/3D/FAULT_SOURCE/branches/decompose_mesh_SCOTCH/Makefile 2011-05-04 01:12:19 UTC (rev 18311)
@@ -0,0 +1,36 @@
+# Makefile
+
+#############################################################
+## Compiler :
+#F90 = ifort # use -warn
+F90 = gfortran # use -Wall
+
+## modify to match your library paths
+SCOTCH_LIBS = -L/SCOTCH_LIB_PATH -lscotch -lscotcherr
+#(SCOTCH library address .Change it in case
+#your SCOTH is install in another address)
+
+#############################################################
+
+LIBS = part_decompose_mesh_SCOTCH.o \
+ decompose_mesh_SCOTCH.o \
+ program_decompose_mesh_SCOTCH.o
+
+# targets
+all: xdecompose_mesh_SCOTCH
+
+xdecompose_mesh_SCOTCH: $(LIBS)
+ ${F90} -o xdecompose_mesh_SCOTCH $(LIBS) $(SCOTCH_LIBS)
+
+
+part_decompose_mesh_SCOTCH.o: part_decompose_mesh_SCOTCH.f90
+ ${F90} -c part_decompose_mesh_SCOTCH.f90
+
+decompose_mesh_SCOTCH.o: decompose_mesh_SCOTCH.f90 part_decompose_mesh_SCOTCH.f90
+ ${F90} -c decompose_mesh_SCOTCH.f90
+
+program_decompose_mesh_SCOTCH.o: program_decompose_mesh_SCOTCH.f90
+ ${F90} -c program_decompose_mesh_SCOTCH.f90
+
+clean:
+ rm -f *.o *.mod a.out xdecompose_mesh_SCOTCH
Added: seismo/3D/FAULT_SOURCE/branches/decompose_mesh_SCOTCH/README
===================================================================
--- seismo/3D/FAULT_SOURCE/branches/decompose_mesh_SCOTCH/README (rev 0)
+++ seismo/3D/FAULT_SOURCE/branches/decompose_mesh_SCOTCH/README 2011-05-04 01:12:19 UTC (rev 18311)
@@ -0,0 +1,38 @@
+----------------------------------------------------------------------
+README
+----------------------------------------------------------------------
+
+
+decompose mesh files:
+
+
+****
+you will need SCOTCH libraries to be installed on your system
+(http://gforge.inria.fr/projects/scotch/)
+
+to compile this executable xdecompose_mesh_SCOTCH for partitioning your mesh files
+****
+
+ 1. create mesh using CUBIT and scripts boundary_definition.py and
+ cubit2specfem3D.py to generate all mesh files
+
+ 2. compile executable "xdecompose_mesh_SCOTCH" in this directory decompose_mesh_SCOTCH/:
+
+ make sure, you have the right location of your SCOTCH library defined
+ in the Makefile (SCOTCH_LIBS), then type:
+
+ > make
+
+ 3. create Database files for the number of partitions/processes you want SPECFEM
+ to run on. These Database files will be needed later for the "xgenerate_databases" executable:
+
+ > ./xdecompose_mesh_SCOTCH nproc input_dir output_dir
+
+ where
+ - nproc is the number of partitions (i.e. parallel processes),
+ - input_dir is the directory containing the mesh files (i.e. "nodes_coord_file","mesh_file",...)
+ - output_dir is the directory to hold the new "proc****_Database" files
+
+ for example, a call could look like:
+
+ > ./xdecompose_mesh_SCOTCH nproc ../CUBIT/MESH/ ../DATABASES_MPI/
Added: seismo/3D/FAULT_SOURCE/branches/decompose_mesh_SCOTCH/compile_all.bash
===================================================================
--- seismo/3D/FAULT_SOURCE/branches/decompose_mesh_SCOTCH/compile_all.bash (rev 0)
+++ seismo/3D/FAULT_SOURCE/branches/decompose_mesh_SCOTCH/compile_all.bash 2011-05-04 01:12:19 UTC (rev 18311)
@@ -0,0 +1,9 @@
+#!/bin/bash
+make clean
+make
+# ./xdecompose_mesh_SCOTCH n input_dir output_dir
+# ./xdecompose_mesh_SCOTCH 4 MESH/OUTPUT_FILES ../DATABASES_MPI/
+npart=`grep nparts constants_decompose_mesh_SCOTCH.h | cut -d = -f 2`
+
+echo "./xdecompose_mesh_SCOTCH $npart $input_dir $ouput_dir "
+./xdecompose_mesh_SCOTCH 'nproc' ../CUBIT/MESH/ ../DATABASES_MPI/
Property changes on: seismo/3D/FAULT_SOURCE/branches/decompose_mesh_SCOTCH/compile_all.bash
___________________________________________________________________
Name: svn:executable
+ *
Added: seismo/3D/FAULT_SOURCE/branches/decompose_mesh_SCOTCH/decompose_mesh_SCOTCH.f90
===================================================================
--- seismo/3D/FAULT_SOURCE/branches/decompose_mesh_SCOTCH/decompose_mesh_SCOTCH.f90 (rev 0)
+++ seismo/3D/FAULT_SOURCE/branches/decompose_mesh_SCOTCH/decompose_mesh_SCOTCH.f90 2011-05-04 01:12:19 UTC (rev 18311)
@@ -0,0 +1,712 @@
+!program pre_meshfem3D
+
+module decompose_mesh_SCOTCH
+
+ use part_decompose_mesh_SCOTCH
+ use fault_scotch
+
+ implicit none
+
+ include './scotchf.h'
+
+! number of partitions
+ integer :: nparts ! e.g. 4 for partitioning for 4 CPUs or 4 processes
+
+! mesh arrays
+ integer(long) :: nspec
+ integer, dimension(:,:), allocatable :: elmnts
+ integer, dimension(:,:), allocatable :: mat
+ integer, dimension(:), allocatable :: part
+
+ integer :: nnodes
+ double precision, dimension(:,:), allocatable :: nodes_coords
+
+ integer, dimension(:), allocatable :: xadj
+ integer, dimension(:), allocatable :: adjncy
+ integer, dimension(:), allocatable :: nnodes_elmnts
+ integer, dimension(:), allocatable :: nodes_elmnts
+ integer, dimension(:), allocatable :: elmnts_load
+
+ integer, dimension(:), pointer :: glob2loc_elmnts
+ integer, dimension(:), pointer :: glob2loc_nodes_nparts
+ integer, dimension(:), pointer :: glob2loc_nodes_parts
+ integer, dimension(:), pointer :: glob2loc_nodes
+
+ integer, dimension(:), pointer :: tab_size_interfaces, tab_interfaces
+ integer, dimension(:), allocatable :: my_interfaces
+ integer, dimension(:), allocatable :: my_nb_interfaces
+ integer :: ninterfaces
+ integer :: my_ninterface
+
+ integer(long) :: nsize ! Max number of elements that contain the same node.
+ integer :: nb_edges
+
+ integer :: ispec, inode
+ integer :: ngnod
+ integer :: max_neighbour ! Real maximum number of neighbours per element
+ integer(long) :: sup_neighbour ! Majoration of the maximum number of neighbours per element
+
+ integer :: ipart, nnodes_loc, nspec_loc
+ integer :: num_elmnt, num_node, num_mat
+
+ ! boundaries
+ integer :: ispec2D
+ integer :: nspec2D_xmin, nspec2D_xmax, nspec2D_ymin, nspec2D_ymax, nspec2D_bottom, nspec2D_top
+ integer, dimension(:), allocatable :: ibelm_xmin, ibelm_xmax, ibelm_ymin, ibelm_ymax, ibelm_bottom, ibelm_top
+ integer, dimension(:,:), allocatable :: nodes_ibelm_xmin, nodes_ibelm_xmax, nodes_ibelm_ymin
+ integer, dimension(:,:), allocatable :: nodes_ibelm_ymax, nodes_ibelm_bottom, nodes_ibelm_top
+
+ ! moho surface (optional)
+ integer :: nspec2D_moho
+ integer, dimension(:), allocatable :: ibelm_moho
+ integer, dimension(:,:), allocatable :: nodes_ibelm_moho
+
+ character(len=256) :: prname
+
+ logical, dimension(:), allocatable :: mask_nodes_elmnts
+ integer, dimension(:), allocatable :: used_nodes_elmnts
+
+ double precision, dimension(SCOTCH_GRAPHDIM) :: scotchgraph
+ double precision, dimension(SCOTCH_STRATDIM) :: scotchstrat
+ character(len=256), parameter :: scotch_strategy='b{job=t,map=t,poli=S,sep=h{pass=30}}'
+ integer :: ierr,idummy
+
+ !pll
+ double precision , dimension(:,:), allocatable :: mat_prop
+ integer :: count_def_mat,count_undef_mat,imat
+ character (len=30), dimension(:,:), allocatable :: undef_mat_prop
+
+! default mesh file directory
+ character(len=256) :: localpath_name
+ character(len=256) :: outputpath_name
+
+ integer :: q_flag,aniso_flag,idomain_id
+ double precision :: vp,vs,rho
+
+ contains
+
+ !----------------------------------------------------------------------------------------------
+ ! reads in mesh files
+ !----------------------------------------------------------------------------------------------
+ subroutine read_mesh_files
+
+ ! sets number of nodes per element (Percy : esize = 8 defined in constants_decom...h)
+ ngnod = esize
+
+ ! reads node coordinates
+ open(unit=98, file=localpath_name(1:len_trim(localpath_name))//'/nodes_coords_file',&
+ status='old', form='formatted', iostat = ierr)
+ if( ierr /= 0 ) then
+ print*,'could not open file:',localpath_name(1:len_trim(localpath_name))//'/nodes_coords_file'
+ stop 'error file open'
+ endif
+ read(98,*) nnodes
+ allocate(nodes_coords(3,nnodes))
+ do inode = 1, nnodes
+ ! format: #id_node #x_coordinate #y_coordinate #z_coordinate
+ read(98,*) num_node, nodes_coords(1,num_node), nodes_coords(2,num_node), nodes_coords(3,num_node)
+ !if(num_node /= inode) stop "ERROR : Invalid nodes_coords file."
+ end do
+ close(98)
+ print*, 'total number of nodes: '
+ print*, ' nnodes = ', nnodes
+
+ ! reads mesh elements indexing
+ !(CUBIT calls this the connectivity, guess in the sense that it connects with the points index in
+ ! the global coordinate file "nodes_coords_file"; it doesn't tell you which point is connected with others)
+ open(unit=98, file=localpath_name(1:len_trim(localpath_name))//'/mesh_file', &
+ status='old', form='formatted')
+ read(98,*) nspec
+ allocate(elmnts(esize,nspec))
+ do ispec = 1, nspec
+ ! format: # element_id #id_node1 ... #id_node8
+
+ ! note: be aware that here we can have different node ordering for a cube element;
+ ! the ordering from Cubit files might not be consistent for multiple volumes, or uneven, unstructured grids
+ !
+ ! guess here it assumes that spectral elements ordering is like first at the bottom of the element, anticlock-wise, i.e.
+ ! point 1 = (0,0,0), point 2 = (0,1,0), point 3 = (1,1,0), point 4 = (1,0,0)
+ ! then top (positive z-direction) of element
+ ! point 5 = (0,0,1), point 6 = (0,1,1), point 7 = (1,1,1), point 8 = (1,0,1)
+
+ !read(98,*) num_elmnt, elmnts(5,num_elmnt), elmnts(1,num_elmnt),elmnts(4,num_elmnt), elmnts(8,num_elmnt), &
+ ! elmnts(6,num_elmnt), elmnts(2,num_elmnt), elmnts(3,num_elmnt), elmnts(7,num_elmnt)
+
+ read(98,*) num_elmnt, elmnts(1,num_elmnt), elmnts(2,num_elmnt),elmnts(3,num_elmnt), elmnts(4,num_elmnt), &
+ elmnts(5,num_elmnt), elmnts(6,num_elmnt), elmnts(7,num_elmnt), elmnts(8,num_elmnt)
+
+ if((num_elmnt > nspec) .or. (num_elmnt < 1) ) stop "ERROR : Invalid mesh file."
+
+
+ !outputs info for each element to see ordering
+ !print*,'ispec: ',ispec
+ !print*,' ',num_elmnt, elmnts(5,num_elmnt), elmnts(1,num_elmnt),elmnts(4,num_elmnt), elmnts(8,num_elmnt), &
+ ! elmnts(6,num_elmnt), elmnts(2,num_elmnt), elmnts(3,num_elmnt), elmnts(7,num_elmnt)
+ !print*,'elem:',num_elmnt
+ !do i=1,8
+ ! print*,' i ',i,'val :',elmnts(i,num_elmnt),&
+ ! nodes_coords(1,elmnts(i,num_elmnt)),nodes_coords(2,elmnts(i,num_elmnt)),nodes_coords(3,elmnts(i,num_elmnt))
+ !enddo
+ !print*
+
+ end do
+ close(98)
+ print*, 'total number of spectral elements:'
+ print*, ' nspec = ', nspec
+
+ ! reads material associations
+ open(unit=98, file=localpath_name(1:len_trim(localpath_name))//'/materials_file', &
+ status='old', form='formatted')
+ allocate(mat(2,nspec))
+ do ispec = 1, nspec
+ ! format: # id_element #flag
+ ! note: be aware that elements may not be sorted in materials_file
+ read(98,*) num_mat,mat(1,num_mat) !mat(1,ispec)!, mat(2,ispec)
+ if((num_mat > nspec) .or. (num_mat < 1) ) stop "ERROR : Invalid mat file."
+ end do
+ close(98)
+
+ ! TODO:
+ ! must be changed, if mat(1,i) < 0 1 == interface , 2 == tomography
+ mat(2,:) = 1
+
+ ! reads material definitions
+ !
+ ! note: format of nummaterial_velocity_file must be
+ !
+ ! #(1)material_domain_id #(2)material_id #(3)rho #(4)vp #(5)vs #(6)Q_flag #(7)anisotropy_flag
+ !
+ ! where
+ ! material_domain_id : 1=acoustic / 2=elastic / 3=poroelastic
+ ! material_id : number of material/volume
+ ! rho : density
+ ! vp : P-velocity
+ ! vs : S-velocity
+ ! Q_flag : 0=no attenuation/1=IATTENUATION_SEDIMENTS_40, 2=..., 13=IATTENUATION_BEDROCK
+ ! anisotropy_flag : 0=no anisotropy/ 1,2,.. check with implementation in aniso_model.f90
+ count_def_mat = 0
+ count_undef_mat = 0
+ open(unit=98, file=localpath_name(1:len_trim(localpath_name))//'/nummaterial_velocity_file',&
+ status='old', form='formatted')
+ ! note: format #material_domain_id #material_id #...
+ read(98,*,iostat=ierr) idummy,num_mat
+ print *,'materials:'
+ ! counts materials (defined/undefined)
+ do while (ierr == 0)
+ print*, ' num_mat = ',num_mat
+ if(num_mat /= -1) then
+ count_def_mat = count_def_mat + 1
+ else
+ count_undef_mat = count_undef_mat + 1
+ end if
+ read(98,*,iostat=ierr) idummy,num_mat
+ end do
+ close(98)
+ print*, ' defined = ',count_def_mat, 'undefined = ',count_undef_mat
+ ! check with material flags
+ if( count_def_mat > 0 .and. maxval(mat(1,:)) > count_def_mat ) then
+ print*,'error material definitions:'
+ print*,' materials associated in materials_file:',maxval(mat(1,:))
+ print*,' bigger than defined materials in nummaterial_velocity_file:',count_def_mat
+ stop 'error materials'
+ endif
+ allocate(mat_prop(6,count_def_mat))
+ allocate(undef_mat_prop(6,count_undef_mat))
+ ! reads in defined material properties
+ open(unit=98, file=localpath_name(1:len_trim(localpath_name))//'/nummaterial_velocity_file', &
+ status='old', form='formatted')
+ do imat=1,count_def_mat
+ ! material definitions
+ !
+ ! format: note that we save the arguments in a slightly different order in mat_prop(:,:)
+ ! #(6) material_domain_id #(0) material_id #(1) rho #(2) vp #(3) vs #(4) Q_flag #(5) anisotropy_flag
+ !
+ read(98,*) idomain_id,num_mat,rho,vp,vs,q_flag,aniso_flag
+ !read(98,*) num_mat, mat_prop(1,num_mat),mat_prop(2,num_mat),&
+ ! mat_prop(3,num_mat),mat_prop(4,num_mat),mat_prop(5,num_mat)
+ mat_prop(1,num_mat) = rho
+ mat_prop(2,num_mat) = vp
+ mat_prop(3,num_mat) = vs
+ mat_prop(4,num_mat) = q_flag
+ mat_prop(5,num_mat) = aniso_flag
+ mat_prop(6,num_mat) = idomain_id
+
+ if(num_mat < 0 .or. num_mat > count_def_mat) stop "ERROR : Invalid nummaterial_velocity_file file."
+
+ !checks attenuation flag with integer range as defined in constants.h like IATTENUATION_SEDIMENTS_40, ....
+ if( int(mat_prop(4,num_mat)) > 13 ) then
+ stop 'wrong attenuation flag in mesh: too large, not supported yet - check with constants.h'
+ endif
+ end do
+ ! reads in undefined material properties
+ do imat=1,count_undef_mat
+ read(98,'(6A30)') undef_mat_prop(6,imat),undef_mat_prop(1,imat),undef_mat_prop(2,imat),&
+ undef_mat_prop(3,imat),undef_mat_prop(4,imat),undef_mat_prop(5,imat)
+ end do
+ close(98)
+
+ ! reads in absorbing boundary files
+ open(unit=98, file=localpath_name(1:len_trim(localpath_name))//'/absorbing_surface_file_xmin', &
+ status='old', form='formatted',iostat=ierr)
+ if( ierr /= 0 ) then
+ nspec2D_xmin = 0
+ else
+ read(98,*) nspec2D_xmin
+ endif
+ allocate(ibelm_xmin(nspec2D_xmin))
+ allocate(nodes_ibelm_xmin(4,nspec2D_xmin))
+ do ispec2D = 1,nspec2D_xmin
+ ! format: #id_(element containing the face) #id_node1_face .. #id_node4_face
+ ! note: ordering for CUBIT seems such that the normal of the face points outward of the element the face belongs to;
+ ! in other words, nodes are in increasing order such that when looking from within the element outwards,
+ ! they are ordered clockwise
+ !
+ ! doesn't necessarily have to start on top-rear, then bottom-rear, bottom-front, and finally top-front i.e.:
+ ! point 1 = (0,1,1), point 2 = (0,1,0), point 3 = (0,0,0), point 4 = (0,0,1)
+ read(98,*) ibelm_xmin(ispec2D), nodes_ibelm_xmin(1,ispec2D), nodes_ibelm_xmin(2,ispec2D), &
+ nodes_ibelm_xmin(3,ispec2D), nodes_ibelm_xmin(4,ispec2D)
+
+ !outputs info for each element for check of ordering
+ !print*,'ispec2d:',ispec2d
+ !print*,' xmin:', ibelm_xmin(ispec2D), nodes_ibelm_xmin(1,ispec2D), nodes_ibelm_xmin(2,ispec2D), &
+ ! nodes_ibelm_xmin(3,ispec2D), nodes_ibelm_xmin(4,ispec2D)
+ !do i=1,4
+ ! print*,'i',i,'val:',ibelm_xmin(ispec2d),nodes_coords(1,nodes_ibelm_xmin(i,ispec2D)), &
+ ! nodes_coords(2,nodes_ibelm_xmin(i,ispec2D)),nodes_coords(3,nodes_ibelm_xmin(i,ispec2D))
+ !enddo
+ !print*
+ end do
+ close(98)
+ print*, 'absorbing boundaries:'
+ print*, ' nspec2D_xmin = ', nspec2D_xmin
+
+ ! reads in absorbing boundary files
+ open(unit=98, file=localpath_name(1:len_trim(localpath_name))//'/absorbing_surface_file_xmax', &
+ status='old', form='formatted',iostat=ierr)
+ if( ierr /= 0 ) then
+ nspec2D_xmax = 0
+ else
+ read(98,*) nspec2D_xmax
+ endif
+ allocate(ibelm_xmax(nspec2D_xmax))
+ allocate(nodes_ibelm_xmax(4,nspec2D_xmax))
+ do ispec2D = 1,nspec2D_xmax
+ ! format: #id_(element containing the face) #id_node1_face .. #id_node4_face
+ read(98,*) ibelm_xmax(ispec2D), nodes_ibelm_xmax(1,ispec2D), nodes_ibelm_xmax(2,ispec2D), &
+ nodes_ibelm_xmax(3,ispec2D), nodes_ibelm_xmax(4,ispec2D)
+ end do
+ close(98)
+ print*, ' nspec2D_xmax = ', nspec2D_xmax
+
+ ! reads in absorbing boundary files
+ open(unit=98, file=localpath_name(1:len_trim(localpath_name))//'/absorbing_surface_file_ymin', &
+ status='old', form='formatted',iostat=ierr)
+ if( ierr /= 0 ) then
+ nspec2D_ymin = 0
+ else
+ read(98,*) nspec2D_ymin
+ endif
+ allocate(ibelm_ymin(nspec2D_ymin))
+ allocate(nodes_ibelm_ymin(4,nspec2D_ymin))
+ do ispec2D = 1,nspec2D_ymin
+ ! format: #id_(element containing the face) #id_node1_face .. #id_node4_face
+ read(98,*) ibelm_ymin(ispec2D), nodes_ibelm_ymin(1,ispec2D), nodes_ibelm_ymin(2,ispec2D), &
+ nodes_ibelm_ymin(3,ispec2D), nodes_ibelm_ymin(4,ispec2D)
+ end do
+ close(98)
+ print*, ' nspec2D_ymin = ', nspec2D_ymin
+
+ ! reads in absorbing boundary files
+ open(unit=98, file=localpath_name(1:len_trim(localpath_name))//'/absorbing_surface_file_ymax', &
+ status='old', form='formatted',iostat=ierr)
+ if( ierr /= 0 ) then
+ nspec2D_ymax = 0
+ else
+ read(98,*) nspec2D_ymax
+ endif
+ allocate(ibelm_ymax(nspec2D_ymax))
+ allocate(nodes_ibelm_ymax(4,nspec2D_ymax))
+ do ispec2D = 1,nspec2D_ymax
+ ! format: #id_(element containing the face) #id_node1_face .. #id_node4_face
+ read(98,*) ibelm_ymax(ispec2D), nodes_ibelm_ymax(1,ispec2D), nodes_ibelm_ymax(2,ispec2D), &
+ nodes_ibelm_ymax(3,ispec2D), nodes_ibelm_ymax(4,ispec2D)
+ end do
+ close(98)
+ print*, ' nspec2D_ymax = ', nspec2D_ymax
+
+ ! reads in absorbing boundary files
+ open(unit=98, file=localpath_name(1:len_trim(localpath_name))//'/absorbing_surface_file_bottom', &
+ status='old', form='formatted',iostat=ierr)
+ if( ierr /= 0 ) then
+ nspec2D_bottom = 0
+ else
+ read(98,*) nspec2D_bottom
+ endif
+ allocate(ibelm_bottom(nspec2D_bottom))
+ allocate(nodes_ibelm_bottom(4,nspec2D_bottom))
+ do ispec2D = 1,nspec2D_bottom
+ ! format: #id_(element containing the face) #id_node1_face .. #id_node4_face
+ read(98,*) ibelm_bottom(ispec2D), nodes_ibelm_bottom(1,ispec2D), nodes_ibelm_bottom(2,ispec2D), &
+ nodes_ibelm_bottom(3,ispec2D), nodes_ibelm_bottom(4,ispec2D)
+ end do
+ close(98)
+ print*, ' nspec2D_bottom = ', nspec2D_bottom
+
+ ! reads in free_surface boundary files
+ open(unit=98, file=localpath_name(1:len_trim(localpath_name))//'/free_surface_file', &
+ status='old', form='formatted',iostat=ierr)
+ if( ierr /= 0 ) then
+ nspec2D_top = 0
+ else
+ read(98,*) nspec2D_top
+ endif
+ allocate(ibelm_top(nspec2D_top))
+ allocate(nodes_ibelm_top(4,nspec2D_top))
+ do ispec2D = 1,nspec2D_top
+ ! format: #id_(element containing the face) #id_node1_face .. #id_node4_face
+ read(98,*) ibelm_top(ispec2D), nodes_ibelm_top(1,ispec2D), nodes_ibelm_top(2,ispec2D), &
+ nodes_ibelm_top(3,ispec2D), nodes_ibelm_top(4,ispec2D)
+ end do
+ close(98)
+ print*, ' nspec2D_top = ', nspec2D_top
+
+ ! reads in moho_surface boundary files (optional)
+ open(unit=98, file=localpath_name(1:len_trim(localpath_name))//'/moho_surface_file', &
+ status='old', form='formatted',iostat=ierr)
+ if( ierr /= 0 ) then
+ nspec2D_moho = 0
+ else
+ read(98,*) nspec2D_moho
+ endif
+ allocate(ibelm_moho(nspec2D_moho))
+ allocate(nodes_ibelm_moho(4,nspec2D_moho))
+ do ispec2D = 1,nspec2D_moho
+ ! format: #id_(element containing the face) #id_node1_face .. #id_node4_face
+ read(98,*) ibelm_moho(ispec2D), nodes_ibelm_moho(1,ispec2D), nodes_ibelm_moho(2,ispec2D), &
+ nodes_ibelm_moho(3,ispec2D), nodes_ibelm_moho(4,ispec2D)
+ end do
+ close(98)
+ if( nspec2D_moho > 0 ) print*, ' nspec2D_moho = ', nspec2D_moho
+
+! Percy & JPA
+
+ ! -----------------------------------
+ ! Reading fault elements
+ ! ----------------------------------
+ call read_fault_files()
+
+
+ end subroutine read_mesh_files
+
+ !----------------------------------------------------------------------------------------------
+ ! checks valence of nodes
+ !----------------------------------------------------------------------------------------------
+
+ subroutine check_valence
+
+ allocate(mask_nodes_elmnts(nnodes))
+ allocate(used_nodes_elmnts(nnodes))
+ mask_nodes_elmnts(:) = .false.
+ used_nodes_elmnts(:) = 0
+ do ispec = 1, nspec
+ do inode = 1, ESIZE
+ mask_nodes_elmnts(elmnts(inode,ispec)) = .true.
+ used_nodes_elmnts(elmnts(inode,ispec)) = used_nodes_elmnts(elmnts(inode,ispec)) + 1
+ enddo
+ enddo
+ print *, 'nodes valence: '
+ print *, ' min = ',minval(used_nodes_elmnts(:)),'max = ', maxval(used_nodes_elmnts(:))
+ do inode = 1, nnodes
+ if (.not. mask_nodes_elmnts(inode)) then
+ stop 'ERROR : nodes not used.'
+ endif
+ enddo
+ nsize = maxval(used_nodes_elmnts(:))
+ sup_neighbour = ngnod * nsize - (ngnod + (ngnod/2 - 1)*nfaces)
+ print*, ' nsize = ',nsize, 'sup_neighbour = ', sup_neighbour
+
+ end subroutine check_valence
+
+ !----------------------------------------------------------------------------------------------
+ ! divides model into partitions using scotch library functions
+ !----------------------------------------------------------------------------------------------
+
+ subroutine scotch_partitioning
+
+ implicit none
+
+ elmnts(:,:) = elmnts(:,:) - 1
+
+ ! determines maximum neighbors based on 1 common node
+ allocate(xadj(1:nspec+1))
+ allocate(adjncy(1:sup_neighbour*nspec))
+ allocate(nnodes_elmnts(1:nnodes))
+ allocate(nodes_elmnts(1:nsize*nnodes))
+ call mesh2dual_ncommonnodes(nspec, nnodes, nsize, sup_neighbour, elmnts, xadj, adjncy, nnodes_elmnts, &
+ nodes_elmnts, max_neighbour, 1)
+ print*, 'mesh2dual: '
+ print*, ' max_neighbour = ',max_neighbour
+
+ nb_edges = xadj(nspec+1)
+
+ ! allocates & initializes partioning of elements
+ allocate(part(1:nspec))
+ part(:) = -1
+
+ ! initializes
+ ! elements load array
+ allocate(elmnts_load(1:nspec))
+
+ ! uniform load
+ elmnts_load(:) = 1
+
+ ! in case of acoustic/elastic simulation, weights elements accordingly
+ call acoustic_elastic_load(elmnts_load,nspec,count_def_mat,mat(1,:),mat_prop)
+
+ ! SCOTCH partitioning
+ call scotchfstratinit (scotchstrat(1), ierr)
+ if (ierr /= 0) then
+ stop 'ERROR : MAIN : Cannot initialize strat'
+ endif
+
+ call scotchfstratgraphmap (scotchstrat(1), trim(scotch_strategy), ierr)
+ if (ierr /= 0) then
+ stop 'ERROR : MAIN : Cannot build strat'
+ endif
+
+ call scotchfgraphinit (scotchgraph (1), ierr)
+ if (ierr /= 0) then
+ stop 'ERROR : MAIN : Cannot initialize graph'
+ endif
+
+ ! fills graph structure : see user manual (scotch_user5.1.pdf, page 72/73)
+ ! arguments: #(1) graph_structure #(2) baseval(either 0/1) #(3) number_of_vertices
+ ! #(4) adjacency_index_array #(5) adjacency_end_index_array (optional)
+ ! #(6) vertex_load_array (optional) #(7) vertex_label_array
+ ! #(7) number_of_arcs #(8) adjacency_array
+ ! #(9) arc_load_array (optional) #(10) ierror
+ call scotchfgraphbuild (scotchgraph (1), 0, nspec, &
+ xadj (1), xadj (1), &
+ elmnts_load (1), xadj (1), &
+ nb_edges, adjncy (1), &
+ adjncy (1), ierr)
+
+ ! w/out element load, but adjacency array
+ !call scotchfgraphbuild (scotchgraph (1), 0, nspec, &
+ ! xadj (1), xadj (1), &
+ ! xadj (1), xadj (1), &
+ ! nb_edges, adjncy (1), &
+ ! adjncy (1), ierr)
+
+
+ if (ierr /= 0) then
+ stop 'ERROR : MAIN : Cannot build graph'
+ endif
+
+ call scotchfgraphcheck (scotchgraph (1), ierr)
+ if (ierr /= 0) then
+ stop 'ERROR : MAIN : Invalid check'
+ endif
+
+ call scotchfgraphpart (scotchgraph (1), nparts, scotchstrat(1),part(1),ierr)
+ if (ierr /= 0) then
+ stop 'ERROR : MAIN : Cannot part graph'
+ endif
+
+ call scotchfgraphexit (scotchgraph (1), ierr)
+ if (ierr /= 0) then
+ stop 'ERROR : MAIN : Cannot destroy graph'
+ endif
+
+ call scotchfstratexit (scotchstrat(1), ierr)
+ if (ierr /= 0) then
+ stop 'ERROR : MAIN : Cannot destroy strat'
+ endif
+
+
+ ! re-partitioning puts acoustic-elastic coupled elements into same partition
+ ! integer :: nfaces_coupled
+ ! integer, dimension(:,:), pointer :: faces_coupled
+ ! call acoustic_elastic_repartitioning (nspec, nnodes, elmnts, &
+ ! count_def_mat, mat(1,:) , mat_prop, &
+ ! sup_neighbour, nsize, &
+ ! nparts, part, nfaces_coupled, faces_coupled)
+
+ ! FAULT-Percy , re-partitioning the fault into same partition .
+ ! integer :: nfaces_coupled
+ ! integer, dimension(:,:), pointer :: faces_coupled
+
+ ! re-partioning puts faults coupled elements into same partition.
+! call faultside1_faultside2_repartitioning (nspec, nnodes, elmnts, &
+! count_def_mat, mat(1,:) , &
+! sup_neighbour, nsize, &
+! nparts, part)
+! FAULT : output part(0) : contains all fault elements
+! : fault_ispec1,fault_ispec2 (fault elements side1 and side2)
+! : fault_iface1,fault_iface2 (fault faces side2 and side2)
+ call fault_collecting_elements(nspec,nnodes,elmnts, &
+ sup_neighbour,esize,nsize,nparts,part)
+
+! re-partitioning puts moho-surface coupled elements into same partition
+ call moho_surface_repartitioning (nspec, nnodes, elmnts, &
+ sup_neighbour, nsize, nparts, part, &
+ nspec2D_moho,ibelm_moho,nodes_ibelm_moho )
+
+
+ ! local number of each element for each partition
+ call Construct_glob2loc_elmnts(nspec, part, glob2loc_elmnts,nparts)
+
+
+ ! local number of each node for each partition
+ call Construct_glob2loc_nodes(nspec, nnodes,nsize, nnodes_elmnts, nodes_elmnts, part, &
+ glob2loc_nodes_nparts, glob2loc_nodes_parts, glob2loc_nodes, nparts)
+
+
+ ! mpi interfaces
+ ! acoustic/elastic boundaries WILL BE SEPARATED into different MPI partitions
+ call Construct_interfaces(nspec, sup_neighbour, part, elmnts, &
+ xadj, adjncy, tab_interfaces, &
+ tab_size_interfaces, ninterfaces, &
+ nparts)
+
+ !or: uncomment if you want acoustic/elastic boundaries NOT to be separated into different MPI partitions
+ !call Construct_interfaces_no_ac_el_sep(nspec, sup_neighbour, part, elmnts, &
+ ! xadj, adjncy, tab_interfaces, &
+ ! tab_size_interfaces, ninterfaces, &
+ ! count_def_mat, mat_prop(3,:), mat(1,:), nparts)
+
+
+ !------------------------------------------------
+ ! close_fault_crack
+ !------------------------------------------------
+
+ call close_faults(nodes_coords,elmnts,nspec,nnodes,esize)
+
+ end subroutine scotch_partitioning
+
+
+ !----------------------------------------------------------------------------------------------
+ ! writes out new Databases files for each partition
+ !----------------------------------------------------------------------------------------------
+
+ subroutine write_mesh_databases
+
+ allocate(my_interfaces(0:ninterfaces-1))
+ allocate(my_nb_interfaces(0:ninterfaces-1))
+
+ ! writes out Database file for each partition
+
+ do ipart = 0, nparts-1
+
+ ! opens output file
+ write(prname, "(i6.6,'_Database')") ipart
+ open(unit=15,file=outputpath_name(1:len_trim(outputpath_name))//'/proc'//prname,&
+ status='unknown', action='write', form='formatted', iostat = ierr)
+ if( ierr /= 0 ) then
+ print*,'error file open:',outputpath_name(1:len_trim(outputpath_name))//'/proc'//prname
+ print*
+ print*,'check if path exists:',outputpath_name(1:len_trim(outputpath_name))
+ stop 'error file open Database'
+ endif
+
+! FAULT : procxxxxx_fault...
+
+ ! opens output file
+ write(prname, "(i6.6,'_Database_fault')") ipart
+ open(unit=16,file=outputpath_name(1:len_trim(outputpath_name))//'/proc'//prname,&
+ status='unknown', action='write', form='formatted', iostat = ierr)
+ if( ierr /= 0 ) then
+ print*,'error file open:',outputpath_name(1:len_trim(outputpath_name))//'/proc'//prname
+ print*
+ print*,'check if path exists:',outputpath_name(1:len_trim(outputpath_name))
+ stop 'error file open Database'
+ endif
+
+
+ ! gets number of nodes
+ call write_glob2loc_nodes_database(15, ipart, nnodes_loc, nodes_coords, &
+ glob2loc_nodes_nparts, glob2loc_nodes_parts, &
+ glob2loc_nodes, nnodes, 1)
+
+! Sticking in space not a both sides Fault nodes.
+
+ ! gets number of spectral elements
+ call write_partition_database(15, ipart, nspec_loc, nspec, elmnts, &
+ glob2loc_elmnts, glob2loc_nodes_nparts, &
+ glob2loc_nodes_parts, glob2loc_nodes, part, mat, ngnod, 1)
+
+! FAULT : get number of fault spectral elements procxxxx_..._fault
+ call write_partion_fault_database(16, ipart, nspec, elmnts, &
+ glob2loc_elmnts,part,1)
+
+
+ ! writes out node coordinate locations
+ write(15,*) nnodes_loc
+
+ call write_glob2loc_nodes_database(15, ipart, nnodes_loc, nodes_coords,&
+ glob2loc_nodes_nparts, glob2loc_nodes_parts, &
+ glob2loc_nodes, nnodes, 2)
+
+ call write_material_properties_database(15,count_def_mat,count_undef_mat, &
+ mat_prop, undef_mat_prop)
+
+ ! writes out spectral element indices
+ write(15,*) nspec_loc
+
+ call write_partition_database(15, ipart, nspec_loc, nspec, elmnts, &
+ glob2loc_elmnts, glob2loc_nodes_nparts, &
+ glob2loc_nodes_parts, glob2loc_nodes, part, mat, ngnod, 2)
+
+
+! FAULT : Writting out procxxxx_..._fault
+ call write_partion_fault_database(16, ipart, nspec, elmnts, &
+ glob2loc_elmnts,part, 2)
+
+
+ ! writes out absorbing/free-surface boundaries
+ call write_boundaries_database(15, ipart, nspec, nspec2D_xmin, nspec2D_xmax, nspec2D_ymin, &
+ nspec2D_ymax, nspec2D_bottom, nspec2D_top, &
+ ibelm_xmin, ibelm_xmax, ibelm_ymin, &
+ ibelm_ymax, ibelm_bottom, ibelm_top, &
+ nodes_ibelm_xmin, nodes_ibelm_xmax, nodes_ibelm_ymin, &
+ nodes_ibelm_ymax, nodes_ibelm_bottom, nodes_ibelm_top, &
+ glob2loc_elmnts, glob2loc_nodes_nparts, &
+ glob2loc_nodes_parts, glob2loc_nodes, part)
+
+ ! gets number of MPI interfaces
+ call Write_interfaces_database(15, tab_interfaces, tab_size_interfaces, ipart, ninterfaces, &
+ my_ninterface, my_interfaces, my_nb_interfaces, &
+ glob2loc_elmnts, glob2loc_nodes_nparts, glob2loc_nodes_parts, &
+ glob2loc_nodes, 1, nparts)
+
+ ! writes out MPI interfaces elements
+ write(15,*) my_ninterface, maxval(my_nb_interfaces)
+
+ call Write_interfaces_database(15, tab_interfaces, tab_size_interfaces, ipart, ninterfaces, &
+ my_ninterface, my_interfaces, my_nb_interfaces, &
+ glob2loc_elmnts, glob2loc_nodes_nparts, glob2loc_nodes_parts, &
+ glob2loc_nodes, 2, nparts)
+
+ ! writes out moho surface (optional)
+ call write_moho_surface_database(15, ipart, nspec, &
+ glob2loc_elmnts, glob2loc_nodes_nparts, &
+ glob2loc_nodes_parts, glob2loc_nodes, part, &
+ nspec2D_moho,ibelm_moho,nodes_ibelm_moho)
+
+ close(15)
+ close(16)
+
+ end do
+ print*, 'partitions: '
+ print*, ' num = ',nparts
+ print*
+ print*, 'Databases files in directory: ',outputpath_name(1:len_trim(outputpath_name))
+ print*, 'finished successfully'
+ print*
+
+ end subroutine write_mesh_databases
+
+!end program pre_meshfem3D
+
+end module decompose_mesh_SCOTCH
+
Added: seismo/3D/FAULT_SOURCE/branches/decompose_mesh_SCOTCH/fault_scotch.f90
===================================================================
--- seismo/3D/FAULT_SOURCE/branches/decompose_mesh_SCOTCH/fault_scotch.f90 (rev 0)
+++ seismo/3D/FAULT_SOURCE/branches/decompose_mesh_SCOTCH/fault_scotch.f90 2011-05-04 01:12:19 UTC (rev 18311)
@@ -0,0 +1,425 @@
+ module fault_scotch
+
+ implicit none
+ type(bc_fault)
+ private
+ integer :: nspec
+ integer, dimension(:), pointer :: ispec1, ispec2, iface1, iface2
+ end type
+
+ type(bc_fault), allocatable, save :: faults(:)
+
+ integer :: nspec_fault_side1, nspec_fault_side2, nspec_fault
+
+ integer , dimension(:),allocatable :: ispec1 ,ispec2, iface1, iface2
+
+ public read_fault_files,fault_collecting_elements
+
+
+ CONTAINS
+!==========================================================================================
+
+ Subroutine read_fault_files
+
+ integer :: nbfaults ,iflt
+
+ open(1,file='../DATA/FAULT/Par_file_faults.in',status='old',action='read',iostat=ier)
+ if (ier==0) then
+ read(1,*) nbfaults
+ allocate(faults(nbfaults))
+ do iflt = 1 , nbfaults
+ call read_single_fault_file(faults(iflt),iflt)
+ enddo
+ else
+ print*, 'Par_file.in has not been found'
+ endif
+ close(1)
+ end subroutine read_fault_files
+
+
+!---------------------------------------------------------------------------------------------------
+
+ Subroutine read_single_fault_file(bcfault,ifault)
+
+! INPUTS :
+ type(bc_fault),intent(inout) :: bcfault
+
+ integer,intent(in) :: ifault
+ character(len=10) :: NTchar
+
+
+ write(NTchar,1) ifault
+ NTchar = adjustl(NTchar)
+1 format(I5)
+
+ open(1,file='../DATA/FAULT/fault_nodes'//NTCHAR//'.dat',status='old',action='read',iostat=ier)
+ if( ier == 0 ) then
+ read(1,*) bcfault%nspec
+ allocate(bcfault%ispec1(bcfault%nspec))
+ allocate(bcfault%ispec2(bcfault%nspec))
+ allocate(bcfault%iface1(bcfault%nspec))
+ allocate(bcfault%iface2(bcfault%nspec))
+ do j=1,bcfault%nspec
+ read(1,*) bcfault%ispec1(j),bcfault%ispec2(j),bcfault%iface1(j),bcfault%iface2(j)
+ enddo
+ else
+ write(6,*) 'fault_nodes.dat does not exit , no fault detected in the domain'
+ endif
+ close(1)
+
+ end Subroutine read_single_fault_file
+
+
+ !--------------------------------------------------
+ ! Repartitioning : two coupled faultside1/side2 elements are transfered to the same partition
+ !--------------------------------------------------
+
+ Subroutine fault_collecting_elements(nelmnts,nnodes,elmnts,sup_neighbour,esize,nsize,nproc,part)
+
+! INPUTS
+ integer(long),intent(in) :: nelmnts,nsize,esize
+ integer(long),intent(in) :: sup_neighbour
+ integer, dimension(0:esize*nelmnts-1),intent(in) :: elmnts
+ integer, intent(in) :: nnodes, nproc
+! OUTPUTS :
+ integer, dimension(0:nelmnts-1),intent(inout) :: part
+! VARIABLES:
+ logical , dimension(nelmnts) :: is_faultside1=.false., &
+ is_faultside2=.false. ! ISPEC1 , ISPEC2.
+ integer :: nbfaults,iflt
+
+ do iflt=1,nbfaults
+ is_faultside1 = .false.
+ is_faultside2 = .false.
+ is_faultside1(faults(iflt)%ispec1) = .true.
+ is_faultside2(faults(iflt)%ispec2) = .true.
+ call faultside1_faultside2_repartitioning (nelmnts, nnodes, elmnts, sup_neighbour, nsize, &
+ nproc, part,is_faultside1,is_faultside2)
+ end do
+
+ end Subroutine fault_collecting_elements
+
+! ---------------------------------------------------------------------------------------------------
+! JPA: we are doing the same steps for side 1 and 2
+! we don' need separate isfault1 isfault2. A single isfault is sufficient
+! To do: reaplce isfault1 and isfault2 by a single isfault
+ Subroutine faultside1_faultside2_repartitioning (nelmnts, nnodes, elmnts, sup_neighbour, nsize, &
+ nproc, part,is_faultside1,is_faultside2)
+
+! INDIVIDUAL FAULT REPARTITION
+
+! part : iproc number of processor partionated. It will altered patching fault elements into the same partion.
+! Part, once is altered , will be input for write_partition_database.
+
+!INPUTS
+ integer(long),intent(in) :: nelmnts
+ integer, intent(in) :: nnodes, nproc
+ integer(long), intent(in) :: sup_neighbour,nsize
+ logical , dimension(nelmnts), intent(in) :: is_faultside1, &
+ is_faultside2 ! ISPEC1 , ISPEC2.
+!OUTPUTS :
+ integer, dimension(0:nelmnts-1),intent(inout) :: part
+
+!LOCAL VARIABLES :
+ integer, dimension(0:esize*nelmnts-1) :: elmnts
+ integer :: nfaces_coupled
+ integer, dimension(0:nelmnts) :: xadj
+ integer, dimension(0:sup_neighbour*nelmnts-1) :: adjncy
+ integer, dimension(0:nnodes-1) :: nnodes_elmnts
+ integer, dimension(0:nsize*nnodes-1) :: nodes_elmnts
+ integer :: max_neighbour
+
+!SHILDING
+ integer :: i,j, iface,iface_coarse,ier,ipart,nproc_null
+ integer :: el, el_adj,el_coarse
+ logical :: is_repartitioned
+ integer, dimension(:), allocatable :: elem_proc_null
+
+
+! downloading processor 0
+ nproc_null = 0
+ do i = 1,nelmnts
+ ! searching for proc = 0 elements
+ if ( part(i) == 0 ) then
+ nproc_null = nproc_null +1
+ end if
+ end do
+
+ print*, 'Elements proc = 0 redistributed in [{nproc}- nproc0] :'
+ print*, nproc_null
+
+ allocate(elem_proc_null(nproc_null))
+
+
+ nproc_null = 0
+ do i = 1,nelmnts
+ ! Filling up proc = 0 elements
+ if ( part(i) == 0 ) then
+ nproc_null = nproc_null +1
+ elem_proc_null(nproc_null) = i
+ end if
+ end do
+
+
+ ! Redistributing proc-0 elements on the rest of processors
+ ipart=0
+ do i = 1, size(elem_proc_null)
+ if ( ipart == nproc ) ipart = 0
+ ipart = ipart +1
+ part(elem_proc_null(i)) = ipart
+ end do
+
+
+! Percy , This is needed to get adjacent element by common face.
+ call mesh2dual_ncommonnodes(nelmnts, nnodes, nsize, sup_neighbour, &
+ elmnts, xadj, adjncy, nnodes_elmnts, &
+ nodes_elmnts, max_neighbour, 4)
+
+
+ ! counts coupled elements
+
+ nfaces_coupled = 0
+ do el = 0, nelmnts-1
+ if ( is_faultside1(el+1) ) then
+ do el_adj = xadj(el), xadj(el+1) - 1
+ if ( is_faultside2((adjncy(el_adj)+1)) ) then
+ nfaces_coupled = nfaces_coupled + 1
+ endif
+ enddo
+ endif
+ enddo
+
+ print*, 'number of fault elements coupled :'
+ print*, nfaces_coupled
+
+ ! coupled elements
+ ! ---------------
+ ! Allocating neighbours with shared fault faces.
+ ! This is done in order to not break MPI interfaces.
+ ! Dirty implementation.
+ ! Fault faces are shield by double coarse neighbours.
+ !
+ ! 1 2
+ ! 1 1 2 2
+ ! 1 1 [1] [2] 2 2
+ ! 1 1 2 2
+ ! 1 2
+
+ ! Allocating elements with double shield coarsing.
+
+
+! Finding ispec_nodes1, ispec_nodes2.
+
+! ispec1_nodes = elmnts(:,ispec1)
+! ispec2_nodes = elmnts(:,ispec2)
+
+
+ do el = 0, nelmnts-1
+ if ( is_faultside1(el+1) ) then
+! side1 = el + 1
+! ispec1_side1 = ispec1(side1)
+
+ do el_adj = xadj(el), xadj(el+1) - 1
+ if (is_faultside2(adjncy(el_adj)+1)) then
+ part(el) = 0
+! side2 = adjncy(el_adj) + 1
+! ispec_side2 = ispec2(side2)
+
+
+
+ do iface = xadj(el), xadj(el+1)-1
+ part(adjncy(iface)) = 0
+ el_coarse = adjncy(iface)
+ do iface_coarse = xadj(el_coarse),xadj(el_coarse+1)-1
+ part(adjncy(iface_coarse)) = 0
+ enddo
+ enddo
+ endif
+ enddo
+ endif
+ enddo
+
+ do el = 0, nelmnts-1
+ if ( is_faultside2(el+1)) then
+ do el_adj = xadj(el), xadj(el+1) - 1
+ if (is_faultside1(adjncy(el_adj)+1)) then
+ part(el) = 0
+ do iface = xadj(el),xadj(el+1)-1
+ part(adjncy(iface)) = 0
+ el_coarse = adjncy(iface)
+ do iface_coarse = xadj(el_coarse),xadj(el_coarse+1)-1
+ part(adjncy(iface_coarse)) = 0
+ enddo
+ enddo
+ endif
+ enddo
+ endif
+ enddo
+ print*, "FAULT SHIELD DOUBLE-COARSE"
+
+
+ end subroutine faultside1_faultside2_repartitioning
+
+! ---------------------------------------------------------------------------------------------------
+
+ subroutine write_fault_partition_database(IIN_database, iproc, nelmnts, elmnts, &
+ glob2loc_elmnts, part, num_phase,esize)
+
+! include './constants_decompose_mesh_SCOTCH.h'
+
+ integer, intent(in) :: IIN_database
+ integer, intent(in) :: num_phase, iproc, esize
+ integer(long), intent(in) :: nelmnts
+ integer, dimension(0:esize*nelmnts-1), intent(in) :: elmnts
+ integer, dimension(0:nelmnts-1), intent(in) :: part,glob2loc_elmnts
+
+ integer :: i,iflt,ispec_fault
+ integer :: nspec_fault_side1,nspec_fault_side2,nspec_fault
+
+ if ( num_phase == 1 ) then
+ ! counts number of fault elements in this partition
+ nspec_fault = 0
+ nspec_fault_side1 = 0
+ nspec_fault_side2 = 0
+ do iflt=1,size(faults)
+ do i = 1,size(faults(iflt)%ispec1)
+ if ( part(faults(iflt)%ispec1(i)-1) == iproc ) then
+ nspec_fault_side1 = nspec_fault_side1 + 1
+
+ endif
+ enddo
+ enddo
+ do iflt=1,size(faults)
+ do i = 1,size(faults(iflt)%ispec2)
+ if ( part(faults(iflt)%ispec2(i)-1) == iproc ) then
+ nspec_fault_side2 = nspec_fault_side2 + 1
+ endif
+ enddo
+ enddo
+ print* , 'ispec1 :'
+ print* , nspec_fault_side1
+ print* , 'ispec2 :'
+ print* , nspec_fault_side2
+
+ if (nspec_fault_side1 == nspec_fault_side2) then
+ nspec_fault = nspec_fault_side1
+ else
+ stop 'Number of fault elements on iproc do not conside'
+ end if
+ allocate(ispec1(nspec_fault))
+ allocate(iface1(nspec_fault))
+ ispec_fault = 0
+ do iflt=1,size(faults)
+ do i = 1,size(faults(iflt)%ispec1)
+ if ( part(faults(iflt)%ispec1(i)-1) == iproc ) then
+ ispec_fault = ispec_fault + 1
+ ispec1(ispec_fault)=glob2loc_elmnts(faults(iflt)%ispec1(i))
+ iface1(ispec_fault)=faults(iflt)%iface1(i)
+
+ endif
+ enddo
+ enddo
+ allocate(ispec2(nspec_fault))
+ allocate(iface2(nspec_fault))
+ ispec_fault = 0
+ do iflt=1,size(faults)
+ do i = 1,size(faults(iflt)%ispec2)
+ if ( part(faults(iflt)%ispec2(i)-1) == iproc ) then
+ ispec_fault = ispec_fault + 1
+ ispec2(ispec_fault)=glob2loc_elmnts(faults(iflt)%ispec2(i))
+ iface2(ispec_fault)=faults(iflt)%iface2(i)
+ endif
+ enddo
+ enddo
+
+ else
+
+ ! Writes ispec1 , ispec2 , iface1 , iface2
+ write(IIN_database,*) nspec_fault
+ do i = 1,nspec_fault
+ write(IIN_database,*) ispec1(i), ispec2(i), iface1(i), iface2(i)
+ enddo
+
+ endif
+
+ end subroutine write_fault_partition_database
+
+! ---------------------------------------------------------------------------------------------------
+!
+
+ subroutine close_faults(nodes_coords,elmnts,nelmnts,nnodes,esize)
+
+ integer ,intent(in) :: nnodes, esize, nelmnts
+ integer, dimension(0:esize*nelmnts-1), intent(in) :: elmnts
+
+ integer, dimension(3,nnodes), intent(in) :: nodes_coords
+
+
+ do iflt=1,size(faults)
+ call close_fault_single(faults(iflt)%ispec1,faults(iflt)%ispec2, &
+ elmnts,nodes_coords,nnodes,esize,nelmnts)
+ enddo
+ end subroutine close_faults
+
+! ---------------------------------------------------------------------------------------------------
+ subroutine close_fault_single(ispec1,ispec2,elmnts,nodes_coords,nnodes,esize,nelmnts)
+
+ integer ,intent(in) :: nnodes, esize, nelmnts
+ integer, dimension(0:esize*nelmnts-1), intent(in) :: elmnts
+ integer , dimension(:), intent(in) :: ispec1,ispec2
+
+ real(kind=CUSTOM_REAL),dimension(3,nnodes),intent(inout) :: nodes_coords
+
+ real(kind=CUSTOM_REAL) :: l,x,y,z,d
+ integer :: iglob1, iglob2, i, j, k1, k2
+
+ logical :: found_it
+
+ l = 1e-3_CUSTOM_REAL
+
+ do i = 1,size(ispec2)
+ do k2=1,esize
+
+ iglob2 = elmnts(k2,ispec2(i))
+ found_it = .false.
+
+ do j = 1,size(ispec1)
+ do k1=1,esize
+
+ iglob1 = elmnts(k1,ispec1(j))
+
+ d = (nodes_coords(1,iglob2)-nodes_coords(1,iglob1))**2 + &
+ (nodes_coords(2,iglob2)-nodes_coords(2,iglob1))**2 + &
+ (nodes_coords(3,iglob2)-nodes_coords(3,iglob1))**2
+ d = sqrt(d)
+
+ if (d <= l) then
+
+ x = (nodes_coords(1,iglob2) + nodes_coords(1,iglob1))/2_CUSTOM_REAL
+ y = (nodes_coords(2,iglob2) + nodes_coords(2,iglob1))/2_CUSTOM_REAL
+ z = (nodes_coords(3,iglob2) + nodes_coords(3,iglob1))/2_CUSTOM_REAL
+
+ nodes_coords(1,iglob2) = x
+ nodes_coords(2,iglob2) = y
+ nodes_coords(3,iglob2) = z
+
+ nodes_coords(1,iglob1) = x
+ nodes_coords(2,iglob1) = y
+ nodes_coords(3,iglob1) = z
+
+ found_it = .true.
+ cycle
+ endif
+
+ enddo
+ if (found_it) cycle
+ enddo
+
+ enddo
+ enddo
+
+ end subroutine close_faults
+
+! ---------------------------------------------------------------------------------------------------
+
+ end module fault_scotch
Added: seismo/3D/FAULT_SOURCE/branches/decompose_mesh_SCOTCH/part_decompose_mesh_SCOTCH.f90
===================================================================
--- seismo/3D/FAULT_SOURCE/branches/decompose_mesh_SCOTCH/part_decompose_mesh_SCOTCH.f90 (rev 0)
+++ seismo/3D/FAULT_SOURCE/branches/decompose_mesh_SCOTCH/part_decompose_mesh_SCOTCH.f90 2011-05-04 01:12:19 UTC (rev 18311)
@@ -0,0 +1,1475 @@
+module part_decompose_mesh_SCOTCH
+
+ implicit none
+
+! Useful kind types
+ integer ,parameter :: short = SELECTED_INT_KIND(4), long = SELECTED_INT_KIND(18)
+
+! Number of nodes per elements.
+ integer, parameter :: ESIZE = 8
+
+! Number of faces per element.
+ integer, parameter :: nfaces = 6
+
+! very large and very small values
+ double precision, parameter :: HUGEVAL = 1.d+30,TINYVAL = 1.d-9
+
+! acoustic-elastic load balancing:
+! assumes that elastic at least ~6 times more expensive than acoustic
+ integer, parameter :: ACOUSTIC_LOAD = 1
+ integer, parameter :: ELASTIC_LOAD = 4
+
+! include './constants_decompose_mesh_SCOTCH.h'
+
+contains
+
+ !-----------------------------------------------
+ ! Creating dual graph (adjacency is defined by 'ncommonnodes' between two elements).
+ !-----------------------------------------------
+ subroutine mesh2dual_ncommonnodes(nelmnts, nnodes, nsize, sup_neighbour, elmnts,&
+ xadj, adjncy, &
+ nnodes_elmnts, nodes_elmnts, &
+ max_neighbour, ncommonnodes)
+
+ integer(long), intent(in) :: nelmnts
+ integer, intent(in) :: nnodes
+ integer(long), intent(in) :: nsize
+ integer(long), intent(in) :: sup_neighbour
+ integer, dimension(0:esize*nelmnts-1), intent(in) :: elmnts
+
+ integer, dimension(0:nelmnts) :: xadj
+ integer, dimension(0:sup_neighbour*nelmnts-1) :: adjncy
+ integer, dimension(0:nnodes-1) :: nnodes_elmnts
+ integer, dimension(0:nsize*nnodes-1) :: nodes_elmnts
+ integer, intent(out) :: max_neighbour
+ integer, intent(in) :: ncommonnodes
+
+ ! local parameters
+ integer :: i, j, k, l, m, nb_edges
+ logical :: is_neighbour
+ integer :: num_node, n
+ integer :: elem_base, elem_target
+ integer :: connectivity
+
+
+ ! initializes
+ xadj(:) = 0
+ adjncy(:) = 0
+ nnodes_elmnts(:) = 0
+ nodes_elmnts(:) = 0
+ nb_edges = 0
+
+ ! list of elements per node
+ do i = 0, esize*nelmnts-1
+ nodes_elmnts(elmnts(i)*nsize+nnodes_elmnts(elmnts(i))) = i/esize
+ nnodes_elmnts(elmnts(i)) = nnodes_elmnts(elmnts(i)) + 1
+ end do
+
+ ! checking which elements are neighbours ('ncommonnodes' criteria)
+ do j = 0, nnodes-1
+ do k = 0, nnodes_elmnts(j)-1
+ do l = k+1, nnodes_elmnts(j)-1
+
+ connectivity = 0
+ elem_base = nodes_elmnts(k+j*nsize)
+ elem_target = nodes_elmnts(l+j*nsize)
+ do n = 1, esize
+ num_node = elmnts(esize*elem_base+n-1)
+ do m = 0, nnodes_elmnts(num_node)-1
+ if ( nodes_elmnts(m+num_node*nsize) == elem_target ) then
+ connectivity = connectivity + 1
+ end if
+ end do
+ end do
+
+ if ( connectivity >= ncommonnodes) then
+
+ is_neighbour = .false.
+
+ do m = 0, xadj(nodes_elmnts(k+j*nsize))
+ if ( .not.is_neighbour ) then
+ if ( adjncy(nodes_elmnts(k+j*nsize)*sup_neighbour+m) == nodes_elmnts(l+j*nsize) ) then
+ is_neighbour = .true.
+
+ end if
+ end if
+ end do
+ if ( .not.is_neighbour ) then
+ adjncy(nodes_elmnts(k+j*nsize)*sup_neighbour+xadj(nodes_elmnts(k+j*nsize))) = nodes_elmnts(l+j*nsize)
+
+ xadj(nodes_elmnts(k+j*nsize)) = xadj(nodes_elmnts(k+j*nsize)) + 1
+ if (xadj(nodes_elmnts(k+j*nsize))>sup_neighbour) stop 'ERROR : too much neighbours per element, modify the mesh.'
+
+ adjncy(nodes_elmnts(l+j*nsize)*sup_neighbour+xadj(nodes_elmnts(l+j*nsize))) = nodes_elmnts(k+j*nsize)
+
+ xadj(nodes_elmnts(l+j*nsize)) = xadj(nodes_elmnts(l+j*nsize)) + 1
+ if (xadj(nodes_elmnts(l+j*nsize))>sup_neighbour) stop 'ERROR : too much neighbours per element, modify the mesh.'
+ end if
+ end if
+ end do
+ end do
+ end do
+
+ max_neighbour = maxval(xadj)
+
+ ! making adjacency arrays compact (to be used for partitioning)
+ do i = 0, nelmnts-1
+ k = xadj(i)
+ xadj(i) = nb_edges
+ do j = 0, k-1
+ adjncy(nb_edges) = adjncy(i*sup_neighbour+j)
+ nb_edges = nb_edges + 1
+ end do
+ end do
+
+ xadj(nelmnts) = nb_edges
+
+
+ end subroutine mesh2dual_ncommonnodes
+
+
+
+ !--------------------------------------------------
+ ! construct local numbering for the elements in each partition
+ !--------------------------------------------------
+ subroutine Construct_glob2loc_elmnts(nelmnts, part, glob2loc_elmnts,nparts)
+
+! include './constants_decompose_mesh_SCOTCH.h'
+
+ integer(long), intent(in) :: nelmnts
+ integer, dimension(0:nelmnts-1), intent(in) :: part
+ integer, dimension(:), pointer :: glob2loc_elmnts
+
+ integer :: num_glob, num_part, nparts
+ integer, dimension(0:nparts-1) :: num_loc
+
+ ! allocates local numbering array
+ allocate(glob2loc_elmnts(0:nelmnts-1))
+
+ ! initializes number of local points per partition
+ do num_part = 0, nparts-1
+ num_loc(num_part) = 0
+ end do
+
+ ! local numbering
+ do num_glob = 0, nelmnts-1
+ ! gets partition
+ num_part = part(num_glob)
+ ! increments local numbering of elements (starting with 0,1,2,...)
+ glob2loc_elmnts(num_glob) = num_loc(num_part)
+ num_loc(num_part) = num_loc(num_part) + 1
+ end do
+
+
+ end subroutine Construct_glob2loc_elmnts
+
+ !--------------------------------------------------
+ ! construct local numbering for the nodes in each partition
+ !--------------------------------------------------
+ subroutine Construct_glob2loc_nodes(nelmnts, nnodes, nsize, nnodes_elmnts, nodes_elmnts, part, &
+ glob2loc_nodes_nparts, glob2loc_nodes_parts, glob2loc_nodes,nparts)
+
+! include './constants_decompose_mesh_SCOTCH.h'
+
+ integer(long), intent(in) :: nelmnts, nsize
+ integer, intent(in) :: nnodes
+ integer, dimension(0:nelmnts-1), intent(in) :: part
+ integer, dimension(0:nnodes-1), intent(in) :: nnodes_elmnts
+ integer, dimension(0:nsize*nnodes-1), intent(in) :: nodes_elmnts
+ integer, dimension(:), pointer :: glob2loc_nodes_nparts
+ integer, dimension(:), pointer :: glob2loc_nodes_parts
+ integer, dimension(:), pointer :: glob2loc_nodes
+
+ integer :: num_node
+ integer :: el
+ integer :: num_part
+ integer :: size_glob2loc_nodes,nparts
+ integer, dimension(0:nparts-1) :: parts_node
+ integer, dimension(0:nparts-1) :: num_parts
+
+ allocate(glob2loc_nodes_nparts(0:nnodes))
+
+ size_glob2loc_nodes = 0
+ parts_node(:) = 0
+
+ do num_node = 0, nnodes-1
+ glob2loc_nodes_nparts(num_node) = size_glob2loc_nodes
+ do el = 0, nnodes_elmnts(num_node)-1
+ parts_node(part(nodes_elmnts(el+nsize*num_node))) = 1
+
+ end do
+
+ do num_part = 0, nparts-1
+ if ( parts_node(num_part) == 1 ) then
+ size_glob2loc_nodes = size_glob2loc_nodes + 1
+ parts_node(num_part) = 0
+
+ end if
+ end do
+
+ end do
+
+ glob2loc_nodes_nparts(nnodes) = size_glob2loc_nodes
+
+ allocate(glob2loc_nodes_parts(0:glob2loc_nodes_nparts(nnodes)-1))
+ allocate(glob2loc_nodes(0:glob2loc_nodes_nparts(nnodes)-1))
+
+ glob2loc_nodes(0) = 0
+
+ parts_node(:) = 0
+ num_parts(:) = 0
+ size_glob2loc_nodes = 0
+
+
+ do num_node = 0, nnodes-1
+ do el = 0, nnodes_elmnts(num_node)-1
+ parts_node(part(nodes_elmnts(el+nsize*num_node))) = 1
+
+ end do
+ do num_part = 0, nparts-1
+
+ if ( parts_node(num_part) == 1 ) then
+ glob2loc_nodes_parts(size_glob2loc_nodes) = num_part
+ glob2loc_nodes(size_glob2loc_nodes) = num_parts(num_part)
+ size_glob2loc_nodes = size_glob2loc_nodes + 1
+ num_parts(num_part) = num_parts(num_part) + 1
+ parts_node(num_part) = 0
+ end if
+
+ end do
+ end do
+
+
+ end subroutine Construct_glob2loc_nodes
+
+
+
+ !--------------------------------------------------
+ ! Construct interfaces between each partitions.
+ ! Two adjacent elements in distinct partitions make an entry in array tab_interfaces :
+ ! 1/ first element, 2/ second element, 3/ number of common nodes, 4/ first node,
+ ! 5/ second node, if relevant.
+
+ ! interface ignores acoustic and elastic elements
+
+ ! Elements with undefined material are considered as elastic elements.
+ !--------------------------------------------------
+ subroutine Construct_interfaces(nelmnts, sup_neighbour, part, elmnts, xadj, adjncy, &
+ tab_interfaces, tab_size_interfaces, ninterfaces, &
+ nparts)
+
+ integer(long), intent(in) :: nelmnts, sup_neighbour
+ integer, dimension(0:nelmnts-1), intent(in) :: part
+ integer, dimension(0:esize*nelmnts-1), intent(in) :: elmnts
+ integer, dimension(0:nelmnts), intent(in) :: xadj
+ integer, dimension(0:sup_neighbour*nelmnts-1), intent(in) :: adjncy
+ integer, dimension(:),pointer :: tab_size_interfaces, tab_interfaces
+ integer, intent(out) :: ninterfaces
+
+ integer, intent(in) :: nparts
+
+ ! local parameters
+ integer :: num_part, num_part_bis, el, el_adj, num_interface, num_edge, ncommon_nodes, &
+ num_node, num_node_bis
+ integer :: i, j
+
+ ! counts number of interfaces between partitions
+ ninterfaces = 0
+ do i = 0, nparts-1
+ do j = i+1, nparts-1
+ ninterfaces = ninterfaces + 1
+ end do
+ end do
+
+ allocate(tab_size_interfaces(0:ninterfaces))
+ tab_size_interfaces(:) = 0
+
+ num_interface = 0
+ num_edge = 0
+
+! determines acoustic/elastic elements based upon given vs velocities
+! and counts same elements for each interface
+ do num_part = 0, nparts-1
+ do num_part_bis = num_part+1, nparts-1
+ do el = 0, nelmnts-1
+ if ( part(el) == num_part ) then
+ ! looks at all neighbor elements
+ do el_adj = xadj(el), xadj(el+1)-1
+ ! adds element if neighbor element lies in next partition
+ if ( part(adjncy(el_adj)) == num_part_bis ) then
+ num_edge = num_edge + 1
+ end if
+
+ end do
+ end if
+ end do
+ ! stores number of elements at interface
+ tab_size_interfaces(num_interface+1) = tab_size_interfaces(num_interface) + num_edge
+ num_edge = 0
+ num_interface = num_interface + 1
+
+ end do
+ end do
+
+
+! stores element indices for elements from above search at each interface
+ num_interface = 0
+ num_edge = 0
+
+ allocate(tab_interfaces(0:(tab_size_interfaces(ninterfaces)*7-1)))
+ tab_interfaces(:) = 0
+
+ do num_part = 0, nparts-1
+ do num_part_bis = num_part+1, nparts-1
+ do el = 0, nelmnts-1
+ if ( part(el) == num_part ) then
+ do el_adj = xadj(el), xadj(el+1)-1
+ ! adds element if in adjacent partition
+ if ( part(adjncy(el_adj)) == num_part_bis ) then
+ tab_interfaces(tab_size_interfaces(num_interface)*7+num_edge*7+0) = el
+ tab_interfaces(tab_size_interfaces(num_interface)*7+num_edge*7+1) = adjncy(el_adj)
+ ncommon_nodes = 0
+ do num_node = 0, esize-1
+ do num_node_bis = 0, esize-1
+ if ( elmnts(el*esize+num_node) == elmnts(adjncy(el_adj)*esize+num_node_bis) ) then
+ tab_interfaces(tab_size_interfaces(num_interface)*7+num_edge*7+3+ncommon_nodes) &
+ = elmnts(el*esize+num_node)
+ ncommon_nodes = ncommon_nodes + 1
+ end if
+ end do
+ end do
+ if ( ncommon_nodes > 0 ) then
+ tab_interfaces(tab_size_interfaces(num_interface)*7+num_edge*7+2) = ncommon_nodes
+ else
+ print *, "Error while building interfaces!", ncommon_nodes
+ end if
+ num_edge = num_edge + 1
+ end if
+ end do
+ end if
+
+ end do
+ num_edge = 0
+ num_interface = num_interface + 1
+ end do
+ end do
+
+ end subroutine Construct_interfaces
+
+
+ !--------------------------------------------------
+ ! Construct interfaces between each partitions.
+ ! Two adjacent elements in distinct partitions make an entry in array tab_interfaces :
+ ! 1/ first element, 2/ second element, 3/ number of common nodes, 4/ first node,
+ ! 5/ second node, if relevant.
+
+ ! No interface between acoustic and elastic elements.
+
+ ! Elements with undefined material are considered as elastic elements.
+ !--------------------------------------------------
+ subroutine Construct_interfaces_no_ac_el_sep(nelmnts, &
+ sup_neighbour, part, elmnts, xadj, adjncy, &
+ tab_interfaces, tab_size_interfaces, ninterfaces, &
+ nb_materials, cs_material, num_material,nparts)
+
+ integer(long), intent(in) :: nelmnts, sup_neighbour
+ integer, dimension(0:nelmnts-1), intent(in) :: part
+ integer, dimension(0:esize*nelmnts-1), intent(in) :: elmnts
+ integer, dimension(0:nelmnts), intent(in) :: xadj
+ integer, dimension(0:sup_neighbour*nelmnts-1), intent(in) :: adjncy
+ integer, dimension(:),pointer :: tab_size_interfaces, tab_interfaces
+ integer, intent(out) :: ninterfaces
+ integer, dimension(1:nelmnts), intent(in) :: num_material
+ ! vs velocities
+ double precision, dimension(1:nb_materials), intent(in) :: cs_material
+
+ integer, intent(in) :: nb_materials,nparts
+
+ ! local parameters
+ integer :: num_part, num_part_bis, el, el_adj, num_interface, num_edge, ncommon_nodes, &
+ num_node, num_node_bis
+ integer :: i, j
+ logical :: is_acoustic_el, is_acoustic_el_adj
+
+ ! counts number of interfaces between partitions
+ ninterfaces = 0
+ do i = 0, nparts-1
+ do j = i+1, nparts-1
+ ninterfaces = ninterfaces + 1
+ end do
+ end do
+
+ allocate(tab_size_interfaces(0:ninterfaces))
+ tab_size_interfaces(:) = 0
+
+ num_interface = 0
+ num_edge = 0
+
+! determines acoustic/elastic elements based upon given vs velocities
+! and counts same elements for each interface
+ do num_part = 0, nparts-1
+ do num_part_bis = num_part+1, nparts-1
+ do el = 0, nelmnts-1
+ if ( part(el) == num_part ) then
+ ! determines whether element is acoustic or not
+ if(num_material(el+1) > 0) then
+ if ( cs_material(num_material(el+1)) < TINYVAL) then
+ is_acoustic_el = .true.
+ else
+ is_acoustic_el = .false.
+ end if
+ else
+ is_acoustic_el = .false.
+ end if
+ ! looks at all neighbor elements
+ do el_adj = xadj(el), xadj(el+1)-1
+ ! determines whether neighbor element is acoustic or not
+ if(num_material(adjncy(el_adj)+1) > 0) then
+ if ( cs_material(num_material(adjncy(el_adj)+1)) < TINYVAL) then
+ is_acoustic_el_adj = .true.
+ else
+ is_acoustic_el_adj = .false.
+ end if
+ else
+ is_acoustic_el_adj = .false.
+ end if
+ ! adds element if neighbor element has same material acoustic/not-acoustic and lies in next partition
+ if ( (part(adjncy(el_adj)) == num_part_bis) .and. (is_acoustic_el .eqv. is_acoustic_el_adj) ) then
+ num_edge = num_edge + 1
+ end if
+ end do
+ end if
+ end do
+ ! stores number of elements at interface
+ tab_size_interfaces(num_interface+1) = tab_size_interfaces(num_interface) + num_edge
+ num_edge = 0
+ num_interface = num_interface + 1
+
+ end do
+ end do
+
+
+! stores element indices for elements from above search at each interface
+ num_interface = 0
+ num_edge = 0
+
+ allocate(tab_interfaces(0:(tab_size_interfaces(ninterfaces)*7-1)))
+ tab_interfaces(:) = 0
+
+ do num_part = 0, nparts-1
+ do num_part_bis = num_part+1, nparts-1
+ do el = 0, nelmnts-1
+ if ( part(el) == num_part ) then
+ if(num_material(el+1) > 0) then
+ if ( cs_material(num_material(el+1)) < TINYVAL) then
+ is_acoustic_el = .true.
+ else
+ is_acoustic_el = .false.
+ end if
+ else
+ is_acoustic_el = .false.
+ end if
+ do el_adj = xadj(el), xadj(el+1)-1
+ if(num_material(adjncy(el_adj)+1) > 0) then
+ if ( cs_material(num_material(adjncy(el_adj)+1)) < TINYVAL) then
+ is_acoustic_el_adj = .true.
+ else
+ is_acoustic_el_adj = .false.
+ end if
+ else
+ is_acoustic_el_adj = .false.
+ end if
+ if ( (part(adjncy(el_adj)) == num_part_bis) .and. (is_acoustic_el .eqv. is_acoustic_el_adj) ) then
+ tab_interfaces(tab_size_interfaces(num_interface)*7+num_edge*7+0) = el
+ tab_interfaces(tab_size_interfaces(num_interface)*7+num_edge*7+1) = adjncy(el_adj)
+ ncommon_nodes = 0
+ do num_node = 0, esize-1
+ do num_node_bis = 0, esize-1
+ if ( elmnts(el*esize+num_node) == elmnts(adjncy(el_adj)*esize+num_node_bis) ) then
+ tab_interfaces(tab_size_interfaces(num_interface)*7+num_edge*7+3+ncommon_nodes) &
+ = elmnts(el*esize+num_node)
+ ncommon_nodes = ncommon_nodes + 1
+ end if
+ end do
+ end do
+ if ( ncommon_nodes > 0 ) then
+ tab_interfaces(tab_size_interfaces(num_interface)*7+num_edge*7+2) = ncommon_nodes
+ else
+ print *, "Error while building interfaces!", ncommon_nodes
+ end if
+ num_edge = num_edge + 1
+ end if
+ end do
+ end if
+
+ end do
+ num_edge = 0
+ num_interface = num_interface + 1
+ end do
+ end do
+
+ end subroutine Construct_interfaces_no_ac_el_sep
+
+
+
+
+
+ !--------------------------------------------------
+ ! Write nodes (their coordinates) pertaining to iproc partition in the corresponding Database
+ !--------------------------------------------------
+ subroutine write_glob2loc_nodes_database(IIN_database, iproc, npgeo, &
+ nodes_coords, glob2loc_nodes_nparts, glob2loc_nodes_parts, &
+ glob2loc_nodes, nnodes, num_phase,ispec1,ispec2)
+
+ integer, intent(in) :: IIN_database
+ integer, intent(in) :: nnodes, iproc, num_phase
+ integer, intent(inout) :: npgeo
+ integer, intent(in),dimension(:) :: ispec1,ispec2
+
+ double precision, dimension(3,nnodes) :: nodes_coords
+ integer, dimension(:), pointer :: glob2loc_nodes_nparts
+ integer, dimension(:), pointer :: glob2loc_nodes_parts
+ integer, dimension(:), pointer :: glob2loc_nodes
+
+ integer :: i, j
+
+ if ( num_phase == 1 ) then
+ ! counts number of points in partition
+ npgeo = 0
+ do i = 0, nnodes-1
+ do j = glob2loc_nodes_nparts(i), glob2loc_nodes_nparts(i+1)-1
+ if ( glob2loc_nodes_parts(j) == iproc ) then
+ npgeo = npgeo + 1
+
+ end if
+
+ end do
+ end do
+ else
+ ! writes out point coordinates
+ 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)
+ end if
+ end do
+ end do
+ end if
+
+ end subroutine Write_glob2loc_nodes_database
+
+
+ !--------------------------------------------------
+ ! Write material properties in the Database
+ !--------------------------------------------------
+ subroutine write_material_properties_database(IIN_database,count_def_mat,count_undef_mat, mat_prop, undef_mat_prop)
+
+ integer, intent(in) :: IIN_database
+ integer, intent(in) :: count_def_mat,count_undef_mat
+ double precision, dimension(6,count_def_mat) :: mat_prop
+ character (len=30), dimension(6,count_undef_mat) :: undef_mat_prop
+ integer :: i
+
+ write(IIN_database,*) count_def_mat,count_undef_mat
+ do i = 1, count_def_mat
+ ! database material definition
+ !
+ ! format: #rho #vp #vs #Q_flag #anisotropy_flag #domain_id
+ !
+ ! (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)
+ 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))
+ end do
+
+ end subroutine write_material_properties_database
+
+
+ !--------------------------------------------------
+ ! Write elements on boundaries (and their four nodes on boundaries) pertaining to iproc partition in the corresponding Database
+ !--------------------------------------------------
+ subroutine write_boundaries_database(IIN_database, iproc, nelmnts, nspec2D_xmin, nspec2D_xmax, &
+ nspec2D_ymin, nspec2D_ymax, nspec2D_bottom, nspec2D_top, &
+ ibelm_xmin, ibelm_xmax, ibelm_ymin, &
+ ibelm_ymax, ibelm_bottom, ibelm_top, &
+ nodes_ibelm_xmin, nodes_ibelm_xmax, nodes_ibelm_ymin, &
+ nodes_ibelm_ymax, nodes_ibelm_bottom, nodes_ibelm_top, &
+ glob2loc_elmnts, glob2loc_nodes_nparts, &
+ glob2loc_nodes_parts, glob2loc_nodes, part )
+
+ integer, intent(in) :: IIN_database
+ integer, intent(in) :: iproc
+ integer(long), intent(in) :: nelmnts
+ integer, intent(in) :: nspec2D_xmin, nspec2D_xmax, nspec2D_ymin, nspec2D_ymax, nspec2D_bottom, nspec2D_top
+ integer, dimension(nspec2D_xmin), intent(in) :: ibelm_xmin
+ integer, dimension(nspec2D_xmax), intent(in) :: ibelm_xmax
+ integer, dimension(nspec2D_ymin), intent(in) :: ibelm_ymin
+ integer, dimension(nspec2D_ymax), intent(in) :: ibelm_ymax
+ integer, dimension(nspec2D_bottom), intent(in) :: ibelm_bottom
+ integer, dimension(nspec2D_top), intent(in) :: ibelm_top
+
+ integer, dimension(4,nspec2D_xmin), intent(in) :: nodes_ibelm_xmin
+ integer, dimension(4,nspec2D_xmax), intent(in) :: nodes_ibelm_xmax
+ integer, dimension(4,nspec2D_ymin), intent(in) :: nodes_ibelm_ymin
+ integer, dimension(4,nspec2D_ymax), intent(in) :: nodes_ibelm_ymax
+ integer, dimension(4,nspec2D_bottom), intent(in) :: nodes_ibelm_bottom
+ integer, dimension(4,nspec2D_top), intent(in) :: nodes_ibelm_top
+ integer, dimension(:), pointer :: glob2loc_elmnts
+ integer, dimension(:), pointer :: glob2loc_nodes_nparts
+ integer, dimension(:), pointer :: glob2loc_nodes_parts
+ integer, dimension(:), pointer :: glob2loc_nodes
+ integer, dimension(1:nelmnts) :: part
+
+ ! local parameters
+ integer :: i,j
+ integer :: loc_node1, loc_node2, loc_node3, loc_node4
+ integer :: loc_nspec2D_xmin,loc_nspec2D_xmax,loc_nspec2D_ymin, &
+ loc_nspec2D_ymax,loc_nspec2D_bottom,loc_nspec2D_top
+
+ ! counts number of elements for boundary at xmin, xmax, ymin, ymax, bottom, top in this partition
+ loc_nspec2D_xmin = 0
+ do i=1,nspec2D_xmin
+ if(part(ibelm_xmin(i)) == iproc) then
+ loc_nspec2D_xmin = loc_nspec2D_xmin + 1
+ end if
+ end do
+ write(IIN_database,*) 1, loc_nspec2D_xmin
+ loc_nspec2D_xmax = 0
+ do i=1,nspec2D_xmax
+ if(part(ibelm_xmax(i)) == iproc) then
+ loc_nspec2D_xmax = loc_nspec2D_xmax + 1
+ end if
+ end do
+ write(IIN_database,*) 2, loc_nspec2D_xmax
+ loc_nspec2D_ymin = 0
+ do i=1,nspec2D_ymin
+ if(part(ibelm_ymin(i)) == iproc) then
+ loc_nspec2D_ymin = loc_nspec2D_ymin + 1
+ end if
+ end do
+ write(IIN_database,*) 3, loc_nspec2D_ymin
+ loc_nspec2D_ymax = 0
+ do i=1,nspec2D_ymax
+ if(part(ibelm_ymax(i)) == iproc) then
+ loc_nspec2D_ymax = loc_nspec2D_ymax + 1
+ end if
+ end do
+ write(IIN_database,*) 4, loc_nspec2D_ymax
+ loc_nspec2D_bottom = 0
+ do i=1,nspec2D_bottom
+ if(part(ibelm_bottom(i)) == iproc) then
+ loc_nspec2D_bottom = loc_nspec2D_bottom + 1
+ end if
+ end do
+ write(IIN_database,*) 5, loc_nspec2D_bottom
+ loc_nspec2D_top = 0
+ do i=1,nspec2D_top
+ if(part(ibelm_top(i)) == iproc) then
+ loc_nspec2D_top = loc_nspec2D_top + 1
+ end if
+ end do
+ write(IIN_database,*) 6, loc_nspec2D_top
+
+ ! outputs element index and element node indices
+ ! note: assumes that element indices in ibelm_* arrays are in the range from 1 to nspec
+ ! (this is assigned by CUBIT, if this changes the following indexing must be changed as well)
+ ! while glob2loc_elmnts(.) is shifted from 0 to nspec-1 thus
+ ! we need to have the arg of glob2loc_elmnts start at 0 ==> glob2loc_nodes(ibelm_** -1 )
+ do i=1,nspec2D_xmin
+ if(part(ibelm_xmin(i)) == iproc) then
+ do j = glob2loc_nodes_nparts(nodes_ibelm_xmin(1,i)-1), glob2loc_nodes_nparts(nodes_ibelm_xmin(1,i))-1
+ if (glob2loc_nodes_parts(j) == iproc ) then
+ loc_node1 = glob2loc_nodes(j)+1
+ end if
+ end do
+ do j = glob2loc_nodes_nparts(nodes_ibelm_xmin(2,i)-1), glob2loc_nodes_nparts(nodes_ibelm_xmin(2,i))-1
+ if (glob2loc_nodes_parts(j) == iproc ) then
+ loc_node2 = glob2loc_nodes(j)+1
+ end if
+ end do
+ do j = glob2loc_nodes_nparts(nodes_ibelm_xmin(3,i)-1), glob2loc_nodes_nparts(nodes_ibelm_xmin(3,i))-1
+ if (glob2loc_nodes_parts(j) == iproc ) then
+ loc_node3 = glob2loc_nodes(j)+1
+ end if
+ end do
+ do j = glob2loc_nodes_nparts(nodes_ibelm_xmin(4,i)-1), glob2loc_nodes_nparts(nodes_ibelm_xmin(4,i))-1
+ if (glob2loc_nodes_parts(j) == iproc ) then
+ 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
+ end if
+ end do
+
+ do i=1,nspec2D_xmax
+ if(part(ibelm_xmax(i)) == iproc) then
+ do j = glob2loc_nodes_nparts(nodes_ibelm_xmax(1,i)-1), glob2loc_nodes_nparts(nodes_ibelm_xmax(1,i))-1
+ if (glob2loc_nodes_parts(j) == iproc ) then
+ loc_node1 = glob2loc_nodes(j)+1
+ end if
+ end do
+ do j = glob2loc_nodes_nparts(nodes_ibelm_xmax(2,i)-1), glob2loc_nodes_nparts(nodes_ibelm_xmax(2,i))-1
+ if (glob2loc_nodes_parts(j) == iproc ) then
+ loc_node2 = glob2loc_nodes(j)+1
+ end if
+ end do
+ do j = glob2loc_nodes_nparts(nodes_ibelm_xmax(3,i)-1), glob2loc_nodes_nparts(nodes_ibelm_xmax(3,i))-1
+ if (glob2loc_nodes_parts(j) == iproc ) then
+ loc_node3 = glob2loc_nodes(j)+1
+ end if
+ end do
+ do j = glob2loc_nodes_nparts(nodes_ibelm_xmax(4,i)-1), glob2loc_nodes_nparts(nodes_ibelm_xmax(4,i))-1
+ if (glob2loc_nodes_parts(j) == iproc ) then
+ 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
+ end if
+ end do
+
+ do i=1,nspec2D_ymin
+ if(part(ibelm_ymin(i)) == iproc) then
+ do j = glob2loc_nodes_nparts(nodes_ibelm_ymin(1,i)-1), glob2loc_nodes_nparts(nodes_ibelm_ymin(1,i))-1
+ if (glob2loc_nodes_parts(j) == iproc ) then
+ loc_node1 = glob2loc_nodes(j)+1
+ end if
+ end do
+ do j = glob2loc_nodes_nparts(nodes_ibelm_ymin(2,i)-1), glob2loc_nodes_nparts(nodes_ibelm_ymin(2,i))-1
+ if (glob2loc_nodes_parts(j) == iproc ) then
+ loc_node2 = glob2loc_nodes(j)+1
+ end if
+ end do
+ do j = glob2loc_nodes_nparts(nodes_ibelm_ymin(3,i)-1), glob2loc_nodes_nparts(nodes_ibelm_ymin(3,i))-1
+ if (glob2loc_nodes_parts(j) == iproc ) then
+ loc_node3 = glob2loc_nodes(j)+1
+ end if
+ end do
+ do j = glob2loc_nodes_nparts(nodes_ibelm_ymin(4,i)-1), glob2loc_nodes_nparts(nodes_ibelm_ymin(4,i))-1
+ if (glob2loc_nodes_parts(j) == iproc ) then
+ 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
+ end if
+ end do
+
+ do i=1,nspec2D_ymax
+ if(part(ibelm_ymax(i)) == iproc) then
+ do j = glob2loc_nodes_nparts(nodes_ibelm_ymax(1,i)-1), glob2loc_nodes_nparts(nodes_ibelm_ymax(1,i))-1
+ if (glob2loc_nodes_parts(j) == iproc ) then
+ loc_node1 = glob2loc_nodes(j)+1
+ end if
+ end do
+ do j = glob2loc_nodes_nparts(nodes_ibelm_ymax(2,i)-1), glob2loc_nodes_nparts(nodes_ibelm_ymax(2,i))-1
+ if (glob2loc_nodes_parts(j) == iproc ) then
+ loc_node2 = glob2loc_nodes(j)+1
+ end if
+ end do
+ do j = glob2loc_nodes_nparts(nodes_ibelm_ymax(3,i)-1), glob2loc_nodes_nparts(nodes_ibelm_ymax(3,i))-1
+ if (glob2loc_nodes_parts(j) == iproc ) then
+ loc_node3 = glob2loc_nodes(j)+1
+ end if
+ end do
+ do j = glob2loc_nodes_nparts(nodes_ibelm_ymax(4,i)-1), glob2loc_nodes_nparts(nodes_ibelm_ymax(4,i))-1
+ if (glob2loc_nodes_parts(j) == iproc ) then
+ 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
+ end if
+ end do
+
+ do i=1,nspec2D_bottom
+ if(part(ibelm_bottom(i)) == iproc) then
+ do j = glob2loc_nodes_nparts(nodes_ibelm_bottom(1,i)-1), glob2loc_nodes_nparts(nodes_ibelm_bottom(1,i))-1
+ if (glob2loc_nodes_parts(j) == iproc ) then
+ loc_node1 = glob2loc_nodes(j)+1
+ end if
+ end do
+ do j = glob2loc_nodes_nparts(nodes_ibelm_bottom(2,i)-1), glob2loc_nodes_nparts(nodes_ibelm_bottom(2,i))-1
+ if (glob2loc_nodes_parts(j) == iproc ) then
+ loc_node2 = glob2loc_nodes(j)+1
+ end if
+ end do
+ do j = glob2loc_nodes_nparts(nodes_ibelm_bottom(3,i)-1), glob2loc_nodes_nparts(nodes_ibelm_bottom(3,i))-1
+ if (glob2loc_nodes_parts(j) == iproc ) then
+ loc_node3 = glob2loc_nodes(j)+1
+ end if
+ end do
+ do j = glob2loc_nodes_nparts(nodes_ibelm_bottom(4,i)-1), glob2loc_nodes_nparts(nodes_ibelm_bottom(4,i))-1
+ if (glob2loc_nodes_parts(j) == iproc ) then
+ 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
+ end if
+ end do
+
+ do i=1,nspec2D_top
+ if(part(ibelm_top(i)) == iproc) then
+ do j = glob2loc_nodes_nparts(nodes_ibelm_top(1,i)-1), glob2loc_nodes_nparts(nodes_ibelm_top(1,i))-1
+ if (glob2loc_nodes_parts(j) == iproc ) then
+ loc_node1 = glob2loc_nodes(j)+1
+ end if
+ end do
+ do j = glob2loc_nodes_nparts(nodes_ibelm_top(2,i)-1), glob2loc_nodes_nparts(nodes_ibelm_top(2,i))-1
+ if (glob2loc_nodes_parts(j) == iproc ) then
+ loc_node2 = glob2loc_nodes(j)+1
+ end if
+ end do
+ do j = glob2loc_nodes_nparts(nodes_ibelm_top(3,i)-1), glob2loc_nodes_nparts(nodes_ibelm_top(3,i))-1
+ if (glob2loc_nodes_parts(j) == iproc ) then
+ loc_node3 = glob2loc_nodes(j)+1
+ end if
+ end do
+ do j = glob2loc_nodes_nparts(nodes_ibelm_top(4,i)-1), glob2loc_nodes_nparts(nodes_ibelm_top(4,i))-1
+ if (glob2loc_nodes_parts(j) == iproc ) then
+ 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
+ end if
+
+ end do
+
+ end subroutine write_boundaries_database
+
+
+ !--------------------------------------------------
+ ! Write elements (their nodes) pertaining to iproc partition in the corresponding Database
+ !--------------------------------------------------
+ subroutine write_partition_database(IIN_database, iproc, nspec, nelmnts, elmnts, &
+ glob2loc_elmnts, glob2loc_nodes_nparts, &
+ glob2loc_nodes_parts, glob2loc_nodes, &
+ part, num_modele, ngnod, num_phase)
+
+! include './constants_decompose_mesh_SCOTCH.h'
+
+ integer, intent(in) :: IIN_database
+ integer, intent(in) :: num_phase, iproc
+ integer(long), intent(in) :: nelmnts
+ integer, intent(inout) :: nspec
+ integer, dimension(0:nelmnts-1) :: part
+ integer, dimension(0:esize*nelmnts-1) :: elmnts
+ integer, dimension(:), pointer :: glob2loc_elmnts
+ integer, dimension(2,nelmnts) :: num_modele
+ integer, dimension(:), pointer :: glob2loc_nodes_nparts
+ integer, dimension(:), pointer :: glob2loc_nodes_parts
+ integer, dimension(:), pointer :: glob2loc_nodes
+ integer, intent(in) :: ngnod
+
+ integer :: i,j,k
+ integer, dimension(0:ngnod-1) :: loc_nodes
+
+ if ( num_phase == 1 ) then
+ ! counts number of spectral elements in this partition
+ nspec = 0
+ do i = 0, nelmnts-1
+ if ( part(i) == iproc ) then
+ nspec = nspec + 1
+ end if
+ end do
+
+ else
+ ! writes out element corner indices
+ do i = 0, nelmnts-1
+ if ( part(i) == iproc ) then
+
+ do j = 0, ngnod-1
+ do k = glob2loc_nodes_nparts(elmnts(i*ngnod+j)), glob2loc_nodes_nparts(elmnts(i*ngnod+j)+1)-1
+
+ if ( glob2loc_nodes_parts(k) == iproc ) then
+ loc_nodes(j) = glob2loc_nodes(k)
+ end if
+ end do
+
+ end do
+
+ ! 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)
+ end if
+ end do
+ end if
+
+
+ end subroutine write_partition_database
+
+
+
+ !--------------------------------------------------
+ ! Write interfaces (element and common nodes) pertaining to iproc partition in the corresponding Database
+ !--------------------------------------------------
+ subroutine write_interfaces_database(IIN_database, tab_interfaces, tab_size_interfaces, iproc, ninterfaces, &
+ my_ninterface, my_interfaces, my_nb_interfaces, glob2loc_elmnts, 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
+ integer, intent(inout) :: my_ninterface
+ integer, dimension(:), pointer :: tab_size_interfaces
+ integer, dimension(:), pointer :: tab_interfaces
+ integer, dimension(0:ninterfaces-1), intent(inout) :: my_interfaces
+ integer, dimension(0:ninterfaces-1), intent(inout) :: my_nb_interfaces
+ integer, dimension(:), pointer :: glob2loc_elmnts
+ integer, dimension(:), pointer :: glob2loc_nodes_nparts
+ integer, dimension(:), pointer :: glob2loc_nodes_parts
+ integer, dimension(:), pointer :: glob2loc_nodes
+
+ integer, dimension(4) :: local_nodes
+ integer :: local_elmnt
+ integer :: num_phase
+
+ integer :: i, j, k, l
+ integer :: num_interface
+
+ integer :: count_faces
+
+ num_interface = 0
+
+ if ( num_phase == 1 ) then
+ ! counts number of interfaces to neighbouring partitions
+ my_interfaces(:) = 0
+ my_nb_interfaces(:) = 0
+
+ ! double loops over all partitions
+ do i = 0, nparts-1
+ do j = i+1, nparts-1
+ ! only counts if specified partition (iproc) appears and interface elements increment
+ if ( (tab_size_interfaces(num_interface) < tab_size_interfaces(num_interface+1)) .and. &
+ (i == iproc .or. j == iproc) ) then
+ ! sets flag
+ my_interfaces(num_interface) = 1
+ ! sets number of elements on interface
+ my_nb_interfaces(num_interface) = tab_size_interfaces(num_interface+1) - tab_size_interfaces(num_interface)
+ end if
+ num_interface = num_interface + 1
+ end do
+ end do
+ my_ninterface = sum(my_interfaces(:))
+
+ else
+ ! writes out MPI interface elements
+ do i = 0, nparts-1
+ 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)
+ else
+ write(IIN_database,*) i, my_nb_interfaces(num_interface)
+ end if
+
+ count_faces = 0
+ do k = tab_size_interfaces(num_interface), tab_size_interfaces(num_interface+1)-1
+ if ( i == iproc ) then
+ local_elmnt = glob2loc_elmnts(tab_interfaces(k*7+0))+1
+ else
+ 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
+ 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, -1, -1
+ case (2)
+ ! edge element
+ 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), -1, -1
+ case (4)
+ ! face element
+ count_faces = count_faces + 1
+ 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
+ do l = glob2loc_nodes_nparts(tab_interfaces(k*7+5)), &
+ glob2loc_nodes_nparts(tab_interfaces(k*7+5)+1)-1
+ if ( glob2loc_nodes_parts(l) == iproc ) then
+ local_nodes(3) = glob2loc_nodes(l)+1
+ end if
+ end do
+ do l = glob2loc_nodes_nparts(tab_interfaces(k*7+6)), &
+ glob2loc_nodes_nparts(tab_interfaces(k*7+6)+1)-1
+ if ( glob2loc_nodes_parts(l) == iproc ) then
+ 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)
+ case default
+ print *, "error in write_interfaces_database!", tab_interfaces(k*7+2), iproc
+ end select
+ end do
+
+ ! outputs infos
+ !print*,' partition MPI interface:',iproc,num_interface
+ !print*,' element faces: ',count_faces
+
+ end if
+
+ num_interface = num_interface + 1
+ end do
+ end do
+
+ end if
+
+ end subroutine write_interfaces_database
+
+ !--------------------------------------------------
+ ! Write elements on surface boundaries (and their four nodes on boundaries)
+ ! pertaining to iproc partition in the corresponding Database
+ !--------------------------------------------------
+ subroutine write_moho_surface_database(IIN_database, iproc, nelmnts, &
+ glob2loc_elmnts, glob2loc_nodes_nparts, &
+ glob2loc_nodes_parts, glob2loc_nodes, part, &
+ nspec2D_moho,ibelm_moho,nodes_ibelm_moho)
+
+ integer, intent(in) :: IIN_database
+ integer, intent(in) :: iproc
+ integer(long), intent(in) :: nelmnts
+
+ integer, dimension(:), pointer :: glob2loc_elmnts
+ integer, dimension(:), pointer :: glob2loc_nodes_nparts
+ integer, dimension(:), pointer :: glob2loc_nodes_parts
+ integer, dimension(:), pointer :: glob2loc_nodes
+ integer, dimension(1:nelmnts) :: part
+
+ integer ,intent(in) :: nspec2D_moho
+ integer ,dimension(nspec2D_moho), intent(in) :: ibelm_moho
+ integer, dimension(4,nspec2D_moho), intent(in) :: nodes_ibelm_moho
+
+ integer :: i,j
+ integer :: loc_node1, loc_node2, loc_node3, loc_node4
+ integer :: loc_nspec2D_moho
+
+ ! counts number of elements for moho surface in this partition
+ ! optional moho
+ loc_nspec2D_moho = 0
+ do i=1,nspec2D_moho
+ if(part(ibelm_moho(i)) == iproc) then
+ loc_nspec2D_moho = loc_nspec2D_moho + 1
+ end if
+ end do
+ ! checks if anything to do
+ if( loc_nspec2D_moho == 0 ) return
+
+ ! format: #surface_id, #number of elements
+ write(IIN_database,*) 7, loc_nspec2D_moho
+
+ ! outputs element index and element node indices
+ ! note: assumes that element indices in ibelm_* arrays are in the range from 1 to nspec
+ ! (this is assigned by CUBIT, if this changes the following indexing must be changed as well)
+ ! while glob2loc_elmnts(.) is shifted from 0 to nspec-1 thus
+ ! we need to have the arg of glob2loc_elmnts start at 0 ==> glob2loc_nodes(ibelm_** -1 )
+
+ ! optional moho
+ do i=1,nspec2D_moho
+ if(part(ibelm_moho(i)) == iproc) then
+ do j = glob2loc_nodes_nparts(nodes_ibelm_moho(1,i)-1), glob2loc_nodes_nparts(nodes_ibelm_moho(1,i))-1
+ if (glob2loc_nodes_parts(j) == iproc ) then
+ loc_node1 = glob2loc_nodes(j)+1
+ end if
+ end do
+ do j = glob2loc_nodes_nparts(nodes_ibelm_moho(2,i)-1), glob2loc_nodes_nparts(nodes_ibelm_moho(2,i))-1
+ if (glob2loc_nodes_parts(j) == iproc ) then
+ loc_node2 = glob2loc_nodes(j)+1
+ end if
+ end do
+ do j = glob2loc_nodes_nparts(nodes_ibelm_moho(3,i)-1), glob2loc_nodes_nparts(nodes_ibelm_moho(3,i))-1
+ if (glob2loc_nodes_parts(j) == iproc ) then
+ loc_node3 = glob2loc_nodes(j)+1
+ end if
+ end do
+ do j = glob2loc_nodes_nparts(nodes_ibelm_moho(4,i)-1), glob2loc_nodes_nparts(nodes_ibelm_moho(4,i))-1
+ if (glob2loc_nodes_parts(j) == iproc ) then
+ 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
+ end if
+
+ end do
+
+ end subroutine write_moho_surface_database
+
+
+
+ !--------------------------------------------------
+ ! loading : sets weights for acoustic/elastic elements to account for different
+ ! expensive calculations in specfem simulations
+ !--------------------------------------------------
+
+ subroutine acoustic_elastic_load (elmnts_load,nelmnts,nb_materials,num_material,mat_prop)
+ !
+ ! note:
+ ! acoustic material = domainID 1 (stored in mat_prop(6,..) )
+ ! elastic material = domainID 2
+ !
+ implicit none
+
+ integer(long),intent(in) :: nelmnts
+ integer, intent(in) :: nb_materials
+
+ ! load weights
+ integer,dimension(1:nelmnts),intent(out) :: elmnts_load
+
+ ! materials
+ integer, dimension(1:nelmnts), intent(in) :: num_material
+ double precision, dimension(6,nb_materials),intent(in) :: mat_prop
+
+ ! local parameters
+ logical, dimension(nb_materials) :: is_acoustic, is_elastic
+ integer :: i,el
+
+ ! sets acoustic/elastic flags for materials
+ is_acoustic(:) = .false.
+ is_elastic(:) = .false.
+ do i = 1, nb_materials
+ ! acoustic material has idomain_id 1
+ if (mat_prop(6,i) == 1 ) then
+ is_acoustic(i) = .true.
+ endif
+ ! elastic material has idomain_id 2
+ if (mat_prop(6,i) == 2 ) then
+ is_elastic(i) = .true.
+ endif
+ enddo
+
+ ! sets weights for elements
+ do el = 0, nelmnts-1
+ ! acoustic element (cheap)
+ if ( is_acoustic(num_material(el+1)) ) then
+ elmnts_load(el+1) = elmnts_load(el+1)*ACOUSTIC_LOAD
+ endif
+ ! elastic element (expensive)
+ if ( is_elastic(num_material(el+1)) ) then
+ elmnts_load(el+1) = elmnts_load(el+1)*ELASTIC_LOAD
+ endif
+ enddo
+
+ end subroutine acoustic_elastic_load
+
+
+ !--------------------------------------------------
+ ! Repartitioning : two coupled acoustic/elastic elements are transfered to the same partition
+ !--------------------------------------------------
+
+ subroutine acoustic_elastic_repartitioning (nelmnts, nnodes, elmnts, &
+ nb_materials, num_material, mat_prop, &
+ sup_neighbour, nsize, &
+ nproc, part, nfaces_coupled, faces_coupled)
+
+ implicit none
+
+ integer(long),intent(in) :: nelmnts
+ integer, intent(in) :: nnodes, nproc, nb_materials
+ integer(long), intent(in) :: sup_neighbour,nsize
+
+ integer, dimension(1:nelmnts), intent(in) :: num_material
+
+ double precision, dimension(6,nb_materials),intent(in) :: mat_prop
+
+ integer, dimension(0:nelmnts-1) :: part
+ integer, dimension(0:esize*nelmnts-1) :: elmnts
+
+ integer, intent(out) :: nfaces_coupled
+ integer, dimension(:,:), pointer :: faces_coupled
+
+
+ logical, dimension(nb_materials) :: is_acoustic, is_elastic
+
+ ! neighbors
+ integer, dimension(:), allocatable :: xadj
+ integer, dimension(:), allocatable :: adjncy
+ integer, dimension(:), allocatable :: nnodes_elmnts
+ integer, dimension(:), allocatable :: nodes_elmnts
+ integer :: max_neighbour
+
+ integer :: i, iface
+ integer :: el, el_adj
+ logical :: is_repartitioned
+
+ ! sets acoustic/elastic flags for materials
+ is_acoustic(:) = .false.
+ is_elastic(:) = .false.
+ do i = 1, nb_materials
+ if (mat_prop(6,i) == 1 ) then
+ is_acoustic(i) = .true.
+ endif
+ if (mat_prop(6,i) == 2 ) then
+ is_elastic(i) = .true.
+ endif
+ enddo
+
+ ! gets neighbors by 4 common nodes (face)
+ allocate(xadj(0:nelmnts))
+ allocate(adjncy(0:sup_neighbour*nelmnts-1))
+ allocate(nnodes_elmnts(0:nnodes-1))
+ allocate(nodes_elmnts(0:nsize*nnodes-1))
+ !call mesh2dual_ncommonnodes(nelmnts, nnodes, elmnts, xadj, adjncy, nnodes_elmnts, nodes_elmnts,4)
+ call mesh2dual_ncommonnodes(nelmnts, nnodes, nsize, sup_neighbour, &
+ elmnts, xadj, adjncy, nnodes_elmnts, &
+ nodes_elmnts, max_neighbour, 4)
+
+ ! counts coupled elements
+ nfaces_coupled = 0
+ do el = 0, nelmnts-1
+ if ( is_acoustic(num_material(el+1)) ) then
+ do el_adj = xadj(el), xadj(el+1) - 1
+ if ( is_elastic(num_material(adjncy(el_adj)+1)) ) then
+ nfaces_coupled = nfaces_coupled + 1
+ endif
+ enddo
+ endif
+ enddo
+
+ ! coupled elements
+ allocate(faces_coupled(2,nfaces_coupled))
+
+ ! stores elements indices
+ nfaces_coupled = 0
+ do el = 0, nelmnts-1
+ if ( is_acoustic(num_material(el+1)) ) then
+ do el_adj = xadj(el), xadj(el+1) - 1
+ if ( is_elastic(num_material(adjncy(el_adj)+1)) ) then
+ nfaces_coupled = nfaces_coupled + 1
+ faces_coupled(1,nfaces_coupled) = el
+ faces_coupled(2,nfaces_coupled) = adjncy(el_adj)
+ endif
+ enddo
+ endif
+ enddo
+
+ ! puts coupled elements into same partition
+ do i = 1, nfaces_coupled*nproc
+ is_repartitioned = .false.
+ do iface = 1, nfaces_coupled
+ if ( part(faces_coupled(1,iface)) /= part(faces_coupled(2,iface)) ) then
+ if ( part(faces_coupled(1,iface)) < part(faces_coupled(2,iface)) ) then
+ part(faces_coupled(2,iface)) = part(faces_coupled(1,iface))
+ else
+ part(faces_coupled(1,iface)) = part(faces_coupled(2,iface))
+ endif
+ is_repartitioned = .true.
+ endif
+ enddo
+ if ( .not. is_repartitioned ) then
+ exit
+ endif
+ enddo
+
+ end subroutine acoustic_elastic_repartitioning
+
+ !--------------------------------------------------
+ ! Repartitioning : two coupled moho surface elements are transfered to the same partition
+ !--------------------------------------------------
+
+ subroutine moho_surface_repartitioning (nelmnts, nnodes, elmnts, &
+ sup_neighbour, nsize, nproc, part, &
+ nspec2D_moho,ibelm_moho,nodes_ibelm_moho)
+
+ implicit none
+
+ ! number of (spectral) elements ( <-> nspec )
+ integer(long),intent(in) :: nelmnts
+
+ ! number of (global) nodes, number or processes
+ integer, intent(in) :: nnodes, nproc
+
+ ! maximum number of neighours and max number of elements-that-contain-the-same-node
+ integer(long), intent(in) :: sup_neighbour,nsize
+
+ ! partition index on each element
+ integer, dimension(0:nelmnts-1) :: part
+
+ ! mesh element indexing
+ ! ( elmnts(esize,nspec) )
+ integer, dimension(0:esize*nelmnts-1) :: elmnts
+
+ ! moho surface
+ integer ,intent(in) :: nspec2D_moho
+ integer ,dimension(nspec2D_moho), intent(in) :: ibelm_moho
+ integer, dimension(4,nspec2D_moho), intent(in) :: nodes_ibelm_moho
+
+ ! local parameters
+ integer :: nfaces_coupled
+ integer, dimension(:,:), pointer :: faces_coupled
+
+ logical, dimension(:),allocatable :: is_moho,node_is_moho
+
+ ! for neighbors
+ integer, dimension(:), allocatable :: xadj
+ integer, dimension(:), allocatable :: adjncy
+ integer, dimension(:), allocatable :: nnodes_elmnts
+ integer, dimension(:), allocatable :: nodes_elmnts
+ integer :: max_neighbour
+
+ integer :: i, j, iface, inode, ispec2D, counter
+ integer :: el, el_adj
+ logical :: is_repartitioned
+
+ ! temporary flag arrays
+ allocate( is_moho(0:nelmnts-1)) ! element ids start from 0
+ allocate( node_is_moho(0:nnodes-1) ) ! node ids start from 0
+ is_moho(:) = .false.
+ node_is_moho(:) = .false.
+
+ ! sets moho flags for known elements
+ do ispec2D = 1, nspec2D_moho
+ ! note: assumes that element indices in ibelm_* arrays are in the range from 1 to nspec
+ el = ibelm_moho(ispec2D) - 1
+ is_moho(el) = .true.
+
+ ! sets node flags
+ do j=1,4
+ ! note: assumes that node indices in nodes_ibelm_* arrays are in the range from 1 to nodes
+ inode = nodes_ibelm_moho(j,ispec2D) - 1
+ node_is_moho(inode) = .true.
+ enddo
+ enddo
+
+ ! checks if element has moho surface
+ do el = 0, nelmnts-1
+ if( is_moho(el) ) cycle
+
+ ! loops over all element corners
+ counter = 0
+ do i=0,esize-1
+ ! note: assumes that node indices in elmnts array are in the range from 0 to nodes-1
+ inode = elmnts(el*esize+i)
+ if( node_is_moho(inode) ) counter = counter + 1
+ enddo
+
+ ! sets flag if it has a surface
+ if( counter == 4 ) is_moho(el) = .true.
+ enddo
+
+ ! statistics output
+ counter = 0
+ do el=0, nelmnts-1
+ if ( is_moho(el) ) counter = counter + 1
+ enddo
+ print*,' moho elements = ',counter
+
+ ! gets neighbors by 4 common nodes (face)
+ allocate(xadj(0:nelmnts)) ! contains number of adjacent elements (neighbours)
+ allocate(adjncy(0:sup_neighbour*nelmnts-1)) ! contains all element id indices of adjacent elements
+ allocate(nnodes_elmnts(0:nnodes-1))
+ allocate(nodes_elmnts(0:nsize*nnodes-1))
+
+ call mesh2dual_ncommonnodes(nelmnts, nnodes, nsize, sup_neighbour, &
+ elmnts, xadj, adjncy, nnodes_elmnts, &
+ nodes_elmnts, max_neighbour, 4)
+
+ ! counts coupled elements
+ nfaces_coupled = 0
+ do el = 0, nelmnts-1
+ if ( is_moho(el) ) then
+ do el_adj = xadj(el), xadj(el+1) - 1
+ ! increments counter if it contains face
+ if( is_moho(adjncy(el_adj)) ) nfaces_coupled = nfaces_coupled + 1
+ enddo
+ endif
+ enddo
+
+ ! coupled elements
+ allocate(faces_coupled(2,nfaces_coupled))
+
+ ! stores elements indices
+ nfaces_coupled = 0
+ do el = 0, nelmnts-1
+ if ( is_moho(el) ) then
+ do el_adj = xadj(el), xadj(el+1) - 1
+ if ( is_moho(adjncy(el_adj)) ) then
+ nfaces_coupled = nfaces_coupled + 1
+ faces_coupled(1,nfaces_coupled) = el
+ faces_coupled(2,nfaces_coupled) = adjncy(el_adj)
+ endif
+ enddo
+ endif
+ enddo
+
+ ! puts coupled elements into same partition
+ do i = 1, nfaces_coupled*nproc
+ is_repartitioned = .false.
+ do iface = 1, nfaces_coupled
+ if ( part(faces_coupled(1,iface)) /= part(faces_coupled(2,iface)) ) then
+ ! coupled moho elements are in different partitions
+ if ( part(faces_coupled(1,iface)) < part(faces_coupled(2,iface)) ) then
+ part(faces_coupled(2,iface)) = part(faces_coupled(1,iface))
+ else
+ part(faces_coupled(1,iface)) = part(faces_coupled(2,iface))
+ endif
+ is_repartitioned = .true.
+ endif
+ enddo
+ if ( .not. is_repartitioned ) then
+ exit
+ endif
+ enddo
+
+ end subroutine moho_surface_repartitioning
+
+end module part_decompose_mesh_SCOTCH
+
Added: seismo/3D/FAULT_SOURCE/branches/decompose_mesh_SCOTCH/program_decompose_mesh_SCOTCH.f90
===================================================================
--- seismo/3D/FAULT_SOURCE/branches/decompose_mesh_SCOTCH/program_decompose_mesh_SCOTCH.f90 (rev 0)
+++ seismo/3D/FAULT_SOURCE/branches/decompose_mesh_SCOTCH/program_decompose_mesh_SCOTCH.f90 2011-05-04 01:12:19 UTC (rev 18311)
@@ -0,0 +1,46 @@
+program pre_meshfem3D
+
+ use decompose_mesh_SCOTCH,only: nparts,localpath_name, outputpath_name,&
+ read_mesh_files, &
+ check_valence, &
+ scotch_partitioning, &
+ write_mesh_databases
+ implicit none
+ integer :: i
+ character(len=256) :: arg(3)
+
+! include './constants_decompose_mesh_SCOTCH.h'
+
+! check usage
+ do i=1,3
+ call getarg(i,arg(i))
+ if (i <= 3 .and. trim(arg(i)) == "") then
+ print *, 'Usage: ./decompose_mesh_SCOTCH nparts input_directory output_directory'
+ print *
+ print *, ' where'
+ print *, ' nparts = number of partitons'
+ print *, ' input_directory = directory containing mesh files mesh_file,nodes_coords_file,..'
+ print *, ' output_directory = directory for output files proc***_Databases'
+ print *
+ stop ' Reenter command line options'
+ endif
+ enddo
+
+ read(arg(1),*) nparts
+ localpath_name = arg(2)
+ outputpath_name = arg(3)
+
+! reads in (CUBIT) mesh files: mesh_file,nodes_coord_file, ...
+ call read_mesh_files()
+
+! checks valence of nodes
+ call check_valence()
+
+! partitions mesh
+ call scotch_partitioning()
+
+! writes out database files
+ call write_mesh_databases()
+
+end program pre_meshfem3D
+
Added: seismo/3D/FAULT_SOURCE/branches/decompose_mesh_SCOTCH/scotchf.h
===================================================================
--- seismo/3D/FAULT_SOURCE/branches/decompose_mesh_SCOTCH/scotchf.h (rev 0)
+++ seismo/3D/FAULT_SOURCE/branches/decompose_mesh_SCOTCH/scotchf.h 2011-05-04 01:12:19 UTC (rev 18311)
@@ -0,0 +1,76 @@
+! ** Copyright 2004,2007 ENSEIRB, INRIA & CNRS
+! **
+! ** This file is part of the Scotch software package for static mapping,
+! ** graph partitioning and sparse matrix ordering.
+! **
+! ** This software is governed by the CeCILL-C license under French law
+! ** and abiding by the rules of distribution of free software. You can
+! ** use, modify and/or redistribute the software under the terms of the
+! ** CeCILL-C license as circulated by CEA, CNRS and INRIA at the following
+! ** URL: "http://www.cecill.info".
+! **
+! ** As a counterpart to the access to the source code and rights to copy,
+! ** modify and redistribute granted by the license, users are provided
+! ** only with a limited warranty and the software's author, the holder of
+! ** the economic rights, and the successive licensors have only limited
+! ** liability.
+! **
+! ** In this respect, the user's attention is drawn to the risks associated
+! ** with loading, using, modifying and/or developing or reproducing the
+! ** software by the user in light of its specific status of free software,
+! ** that may mean that it is complicated to manipulate, and that also
+! ** therefore means that it is reserved for developers and experienced
+! ** professionals having in-depth computer knowledge. Users are therefore
+! ** encouraged to load and test the software's suitability as regards
+! ** their requirements in conditions enabling the security of their
+! ** systems and/or data to be ensured and, more generally, to use and
+! ** operate it in the same conditions as regards security.
+! **
+! ** The fact that you are presently reading this means that you have had
+! ** knowledge of the CeCILL-C license and that you accept its terms.
+! **
+! ************************************************************
+! ** **
+! ** NAME : ptscotchf.h **
+! ** **
+! ** AUTHOR : Francois PELLEGRINI **
+! ** **
+! ** FUNCTION : FORTRAN declaration file for the **
+! ** LibScotch static mapping and sparse **
+! ** matrix block ordering sequential **
+! ** library. **
+! ** **
+! ** DATES : # Version 3.4 : from : 04 feb 2000 **
+! ** to 22 oct 2001 **
+! ** # Version 4.0 : from : 16 jan 2004 **
+! ** to 16 jan 2004 **
+! ** # Version 5.0 : from : 26 apr 2006 **
+! ** to 26 apr 2006 **
+! ** **
+! ************************************************************
+
+! ** Size definitions for the SCOTCH opaque
+! ** structures. These structures must be
+! ** allocated as arrays of DOUBLEPRECISION
+! ** values for proper padding. The dummy
+! ** sizes are computed at compile-time by
+! ** program "dummysizes".
+
+ INTEGER SCOTCH_ARCHDIM
+ INTEGER SCOTCH_DGRAPHDIM
+ INTEGER SCOTCH_DORDERDIM
+ INTEGER SCOTCH_GEOMDIM
+ INTEGER SCOTCH_GRAPHDIM
+ INTEGER SCOTCH_MAPDIM
+ INTEGER SCOTCH_MESHDIM
+ INTEGER SCOTCH_ORDERDIM
+ INTEGER SCOTCH_STRATDIM
+ PARAMETER (SCOTCH_ARCHDIM = 5)
+ PARAMETER (SCOTCH_DGRAPHDIM = 1)
+ PARAMETER (SCOTCH_DORDERDIM = 1)
+ PARAMETER (SCOTCH_GEOMDIM = 2)
+ PARAMETER (SCOTCH_GRAPHDIM = 12)
+ PARAMETER (SCOTCH_MAPDIM = 13)
+ PARAMETER (SCOTCH_MESHDIM = 15)
+ PARAMETER (SCOTCH_ORDERDIM = 12)
+ PARAMETER (SCOTCH_STRATDIM = 1)
More information about the CIG-COMMITS
mailing list