[cig-commits] r14302 - in seismo/3D/SPECFEM3D_SESAME/trunk: . UTILS/external_mesh UTILS/external_mesh/decimate_mesh UTILS/external_mesh/decompose_mesh UTILS/external_mesh/files_needed_asteroid analyze_CUBIT_Abaqus_mesh
dkomati1 at geodynamics.org
dkomati1 at geodynamics.org
Thu Mar 12 10:42:52 PDT 2009
Author: dkomati1
Date: 2009-03-12 10:42:51 -0700 (Thu, 12 Mar 2009)
New Revision: 14302
Added:
seismo/3D/SPECFEM3D_SESAME/trunk/UTILS/external_mesh/TODO_list_external_mesh
seismo/3D/SPECFEM3D_SESAME/trunk/UTILS/external_mesh/decimate_mesh/
seismo/3D/SPECFEM3D_SESAME/trunk/UTILS/external_mesh/decimate_mesh/decimate_mesh.f90
seismo/3D/SPECFEM3D_SESAME/trunk/UTILS/external_mesh/decompose_mesh/
seismo/3D/SPECFEM3D_SESAME/trunk/UTILS/external_mesh/decompose_mesh/compile_all.csh
seismo/3D/SPECFEM3D_SESAME/trunk/UTILS/external_mesh/decompose_mesh/constants_decompose_mesh.h
seismo/3D/SPECFEM3D_SESAME/trunk/UTILS/external_mesh/decompose_mesh/decompose_mesh.f90
seismo/3D/SPECFEM3D_SESAME/trunk/UTILS/external_mesh/decompose_mesh/part_decompose_mesh.f90
seismo/3D/SPECFEM3D_SESAME/trunk/UTILS/external_mesh/files_needed_asteroid/
seismo/3D/SPECFEM3D_SESAME/trunk/UTILS/external_mesh/readme_external_mesh
seismo/3D/SPECFEM3D_SESAME/trunk/analyze_CUBIT_Abaqus_mesh/
Removed:
seismo/3D/SPECFEM3D_SESAME/trunk/UTILS/external_mesh/TODO_list
seismo/3D/SPECFEM3D_SESAME/trunk/UTILS/external_mesh/analyze_CUBIT_Abaqus_mesh/
seismo/3D/SPECFEM3D_SESAME/trunk/UTILS/external_mesh/decimate_mesh/sub.f90
seismo/3D/SPECFEM3D_SESAME/trunk/UTILS/external_mesh/decompose_mesh/compil_all.sh
seismo/3D/SPECFEM3D_SESAME/trunk/UTILS/external_mesh/decompose_mesh/constants_pre_meshfem3D.h
seismo/3D/SPECFEM3D_SESAME/trunk/UTILS/external_mesh/decompose_mesh/part_pre_meshfem3D.f90
seismo/3D/SPECFEM3D_SESAME/trunk/UTILS/external_mesh/decompose_mesh/pre_meshfem3D.f90
seismo/3D/SPECFEM3D_SESAME/trunk/UTILS/external_mesh/files_needed_for_asteroid_simulation/
seismo/3D/SPECFEM3D_SESAME/trunk/UTILS/external_mesh/model_asteroid_subdivide/
seismo/3D/SPECFEM3D_SESAME/trunk/UTILS/external_mesh/pre_meshfem3D_asteroid_subdivide/
seismo/3D/SPECFEM3D_SESAME/trunk/UTILS/external_mesh/readme
Log:
renamed some directories and files
Deleted: seismo/3D/SPECFEM3D_SESAME/trunk/UTILS/external_mesh/TODO_list
===================================================================
--- seismo/3D/SPECFEM3D_SESAME/trunk/UTILS/external_mesh/TODO_list 2009-03-12 17:37:00 UTC (rev 14301)
+++ seismo/3D/SPECFEM3D_SESAME/trunk/UTILS/external_mesh/TODO_list 2009-03-12 17:42:51 UTC (rev 14302)
@@ -1,29 +0,0 @@
-
-DESIGN :
-
- - The subdivision of the mesh should be done inside meshfem3D (there is actually no point in doing it prior to pre_meshfem). The subdivision itself is already implemented, but we also have to subdivide the communications and any other similar data that concerns elements (absorbing boundaries, PML,...).
-
- - Nearly all constants declared in constants.h concerning external meshes (look for string 'nlegoff') should be read/computed from a kind of Par_file. The same goes for the declaration of materials, done with a function in meshfem3D.f90 (and also right now in specfem3D.f90 in case of associating materials after pre_meshfem but we should get rid of it). There are also some variables declared as parameters used in create_movie_AVS_DX.f90.
-
- - Materials should be associated to elements before pre_meshfem, but it is not always possible (that was the case for Celine Blitz's asteroid mesh with regolith because of limitations in the CUBIT software with huge meshes). An example of how we can isolate elements in contact with the envelope and associate them with another material is commented in specfem3D.f90 (look for string 'NL NL REGOLITH') tough it should in fact be inside meshfem3D, not specfem3D. Identifying the envelope is based on the NGLL points, so there is no problem with tranfering it to meshfem3D.
-
- - Calculations on elements in specfem3D for regular meshes should be merged inside compute_forces (the basic calculation is okay, but not attenuation, absorbing boundaries, ...) as there is no reason to differenciate those.
-
-MISC :
-
- - The following default options in flags.guess -qlanglvl=2003pure (for xlf) and -std=f2003 (for gfortran) causes some problems with intrinsics. There might be a better way than to simply remove these options.
-
- - Reading the file DATA/STATIONS_FILTERED with pgf90 returned an error, because it considered that the lines were truncated, thus failed to read some data even when DATA/STATIONS_FILTERED was also generated by pgf90. Go figures. There might be a better workaround than simply removing writing of blanks ' ' (see rev).
-
-Add-ONS :
-
- - Receivers distance to source along the surface computed for an asteroid can be done by computing the shortest path in a graph. No problem in serial, doing so in parallel might be tricky.
-
- - Partitioning should be done in meshfem3D, using a parallel partitioner such as ParMETIS and/or SCOTCH (PT-SCOTCH). Reading the (distributed or not) mesh, building the graph, distributing, and partitioning is not a problem, but redistributing the graph/mesh according to the partitioning obtained might well be.
-
- - Adding absorbing boundaries, attenuation, PML, moment sources, ... for external meshes.
-
- - Differentiate between elements in contact with the envelope, and those that are in contact with a fracture : source, receivers, ... anything based on the envelope detection will have to take into account those fractures. Really tricky without prior information on the mesh.
-
- - Settle the issue of normal recording for the receivers. The normal is unique, but the plane can be described using two vectors at random (while conserving an orthonormal reference), so we have no reason to choose a pair of vectors compared to the others.i
-
Copied: seismo/3D/SPECFEM3D_SESAME/trunk/UTILS/external_mesh/TODO_list_external_mesh (from rev 14301, seismo/3D/SPECFEM3D_SESAME/trunk/UTILS/external_mesh/TODO_list)
===================================================================
--- seismo/3D/SPECFEM3D_SESAME/trunk/UTILS/external_mesh/TODO_list_external_mesh (rev 0)
+++ seismo/3D/SPECFEM3D_SESAME/trunk/UTILS/external_mesh/TODO_list_external_mesh 2009-03-12 17:42:51 UTC (rev 14302)
@@ -0,0 +1,29 @@
+
+DESIGN :
+
+ - The subdivision of the mesh should be done inside meshfem3D (there is actually no point in doing it prior to pre_meshfem). The subdivision itself is already implemented, but we also have to subdivide the communications and any other similar data that concerns elements (absorbing boundaries, PML,...).
+
+ - Nearly all constants declared in constants.h concerning external meshes (look for string 'nlegoff') should be read/computed from a kind of Par_file. The same goes for the declaration of materials, done with a function in meshfem3D.f90 (and also right now in specfem3D.f90 in case of associating materials after pre_meshfem but we should get rid of it). There are also some variables declared as parameters used in create_movie_AVS_DX.f90.
+
+ - Materials should be associated to elements before pre_meshfem, but it is not always possible (that was the case for Celine Blitz's asteroid mesh with regolith because of limitations in the CUBIT software with huge meshes). An example of how we can isolate elements in contact with the envelope and associate them with another material is commented in specfem3D.f90 (look for string 'NL NL REGOLITH') tough it should in fact be inside meshfem3D, not specfem3D. Identifying the envelope is based on the NGLL points, so there is no problem with tranfering it to meshfem3D.
+
+ - Calculations on elements in specfem3D for regular meshes should be merged inside compute_forces (the basic calculation is okay, but not attenuation, absorbing boundaries, ...) as there is no reason to differenciate those.
+
+MISC :
+
+ - The following default options in flags.guess -qlanglvl=2003pure (for xlf) and -std=f2003 (for gfortran) causes some problems with intrinsics. There might be a better way than to simply remove these options.
+
+ - Reading the file DATA/STATIONS_FILTERED with pgf90 returned an error, because it considered that the lines were truncated, thus failed to read some data even when DATA/STATIONS_FILTERED was also generated by pgf90. Go figures. There might be a better workaround than simply removing writing of blanks ' ' (see rev).
+
+Add-ONS :
+
+ - Receivers distance to source along the surface computed for an asteroid can be done by computing the shortest path in a graph. No problem in serial, doing so in parallel might be tricky.
+
+ - Partitioning should be done in meshfem3D, using a parallel partitioner such as ParMETIS and/or SCOTCH (PT-SCOTCH). Reading the (distributed or not) mesh, building the graph, distributing, and partitioning is not a problem, but redistributing the graph/mesh according to the partitioning obtained might well be.
+
+ - Adding absorbing boundaries, attenuation, PML, moment sources, ... for external meshes.
+
+ - Differentiate between elements in contact with the envelope, and those that are in contact with a fracture : source, receivers, ... anything based on the envelope detection will have to take into account those fractures. Really tricky without prior information on the mesh.
+
+ - Settle the issue of normal recording for the receivers. The normal is unique, but the plane can be described using two vectors at random (while conserving an orthonormal reference), so we have no reason to choose a pair of vectors compared to the others.i
+
Property changes on: seismo/3D/SPECFEM3D_SESAME/trunk/UTILS/external_mesh/TODO_list_external_mesh
___________________________________________________________________
Name: svn:mergeinfo
+
Copied: seismo/3D/SPECFEM3D_SESAME/trunk/UTILS/external_mesh/decimate_mesh (from rev 14300, seismo/3D/SPECFEM3D_SESAME/trunk/UTILS/external_mesh/model_asteroid_subdivide)
Property changes on: seismo/3D/SPECFEM3D_SESAME/trunk/UTILS/external_mesh/decimate_mesh
___________________________________________________________________
Name: svn:mergeinfo
+
Copied: seismo/3D/SPECFEM3D_SESAME/trunk/UTILS/external_mesh/decimate_mesh/decimate_mesh.f90 (from rev 14301, seismo/3D/SPECFEM3D_SESAME/trunk/UTILS/external_mesh/model_asteroid_subdivide/decimate_mesh.f90)
===================================================================
--- seismo/3D/SPECFEM3D_SESAME/trunk/UTILS/external_mesh/decimate_mesh/decimate_mesh.f90 (rev 0)
+++ seismo/3D/SPECFEM3D_SESAME/trunk/UTILS/external_mesh/decimate_mesh/decimate_mesh.f90 2009-03-12 17:42:51 UTC (rev 14302)
@@ -0,0 +1,594 @@
+
+ program decimate_mesh
+
+! mesh decimation, by Nicolas Le Goff, 2008
+
+! cuts each hexahedron of a given mesh in 8 hexahedra (2 * 2 * 2)
+! recursively in order to increase mesh density by a factor of 2
+
+ implicit none
+
+ include './constants.h'
+
+ integer :: nelmnts_ext_mesh
+ integer, dimension(:,:), allocatable :: elmnts_ext_mesh
+ integer, dimension(:,:), allocatable :: elmnts_ext_mesh_sub
+ integer, dimension(:), allocatable :: mat_ext_mesh
+ integer, dimension(:), allocatable :: mat_ext_mesh_sub
+
+ integer :: nnodes_ext_mesh
+ real, dimension(:,:), allocatable :: nodes_coords_ext_mesh
+ real, dimension(:,:), allocatable :: nodes_coords_ext_mesh_sub
+
+ real, dimension(NDIM,NSUB+1,NSUB+1,NSUB+1) :: temporary_nodes
+ integer, dimension(NSUB+1,NSUB+1,NSUB+1) :: temporary_nodes_lookup
+
+ integer, dimension(:), allocatable :: xadj
+ integer, dimension(:), allocatable :: adjncy
+ integer, dimension(:), allocatable :: nnodes_elmnts
+ integer, dimension(:), allocatable :: nodes_elmnts
+
+ integer :: ispec, inode, ispec_neighbours, ispec_neighbours_sub
+ integer :: nnodes_ext_mesh_sub
+ integer :: i, j, k
+ integer :: ix, iy, iz
+ integer :: idim
+
+ real :: xtol
+
+ real :: xminval,yminval,xmaxval,ymaxval,xtypdist,zminval,zmaxval
+
+ open(unit=98, file='./mesh', status='old', form='formatted')
+ read(98,*) nelmnts_ext_mesh
+ allocate(elmnts_ext_mesh(ESIZE,nelmnts_ext_mesh))
+ do ispec = 1, nelmnts_ext_mesh
+ read(98,*) elmnts_ext_mesh(1,ispec), elmnts_ext_mesh(2,ispec), elmnts_ext_mesh(3,ispec), elmnts_ext_mesh(4,ispec), &
+ elmnts_ext_mesh(5,ispec), elmnts_ext_mesh(6,ispec), elmnts_ext_mesh(7,ispec), elmnts_ext_mesh(8,ispec)
+ end do
+ close(98)
+
+ open(unit=98, file='./mat', status='old', form='formatted')
+ allocate(mat_ext_mesh(nelmnts_ext_mesh))
+ do ispec = 1, nelmnts_ext_mesh
+ read(98,*) mat_ext_mesh(ispec)
+ end do
+ close(98)
+
+ open(unit=98, file='./nodes_coords', status='old', form='formatted')
+ read(98,*) nnodes_ext_mesh
+ allocate(nodes_coords_ext_mesh(NDIM,nnodes_ext_mesh))
+ do inode = 1, nnodes_ext_mesh
+ read(98,*) nodes_coords_ext_mesh(1,inode), nodes_coords_ext_mesh(2,inode),nodes_coords_ext_mesh(3,inode)
+ end do
+ close(98)
+
+
+! check that there really are 8 nodes per element.
+ do ispec = 1, nelmnts_ext_mesh
+ do inode = 1, ESIZE
+ do ix = inode+1, ESIZE
+ if (elmnts_ext_mesh(inode,ispec) == elmnts_ext_mesh(ix,ispec)) then
+ stop 'ERRORERROR'
+ endif
+ enddo
+
+ enddo
+enddo
+
+! set up local geometric tolerances
+ xtypdist=+HUGEVAL
+
+ do ispec=1,nelmnts_ext_mesh
+
+ xminval=+HUGEVAL
+ yminval=+HUGEVAL
+ zminval=+HUGEVAL
+ xmaxval=-HUGEVAL
+ ymaxval=-HUGEVAL
+ zmaxval=-HUGEVAL
+
+ do inode = 1, 8
+ xmaxval=max(nodes_coords_ext_mesh(1,elmnts_ext_mesh(inode,ispec)),xmaxval)
+ xminval=min(nodes_coords_ext_mesh(1,elmnts_ext_mesh(inode,ispec)),xminval)
+ ymaxval=max(nodes_coords_ext_mesh(2,elmnts_ext_mesh(inode,ispec)),ymaxval)
+ yminval=min(nodes_coords_ext_mesh(2,elmnts_ext_mesh(inode,ispec)),yminval)
+ zmaxval=max(nodes_coords_ext_mesh(3,elmnts_ext_mesh(inode,ispec)),zmaxval)
+ zminval=min(nodes_coords_ext_mesh(3,elmnts_ext_mesh(inode,ispec)),zminval)
+ enddo
+
+! compute the minimum typical "size" of an element in the mesh
+ xtypdist = min(xtypdist,xmaxval-xminval)
+ xtypdist = min(xtypdist,ymaxval-yminval)
+ xtypdist = min(xtypdist,zmaxval-zminval)
+ !xtypdist = min(xtypdist,sqrt((xmaxval-xminval)**2 + (ymaxval-yminval)**2 + (zmaxval-zminval)**2))
+
+ enddo
+
+! define a tolerance, small with respect to the minimum size
+ xtol=smallval_tol*xtypdist*1.d7
+
+ print *, 'xtypdist' , xtypdist
+ print *, 'facteur de tolerance XTOL = ', xtol
+
+ print *, 'xmin', minval(nodes_coords_ext_mesh(1,:))
+ print *, 'xmax', maxval(nodes_coords_ext_mesh(1,:))
+ print *, 'ymin', minval(nodes_coords_ext_mesh(2,:))
+ print *, 'ymax', maxval(nodes_coords_ext_mesh(2,:))
+ print *, 'zmin', minval(nodes_coords_ext_mesh(3,:))
+ print *, 'zmax', maxval(nodes_coords_ext_mesh(3,:))
+
+
+
+! we build the graph
+ elmnts_ext_mesh(:,:) = elmnts_ext_mesh(:,:) - 1
+
+ allocate(xadj(1:nelmnts_ext_mesh+1))
+ allocate(adjncy(1:MAX_NEIGHBOURS*nelmnts_ext_mesh))
+ allocate(nnodes_elmnts(1:nnodes_ext_mesh))
+ allocate(nodes_elmnts(1:NSIZE*nnodes_ext_mesh))
+
+ call mesh2dual_ncommonnodes(nelmnts_ext_mesh, nnodes_ext_mesh, elmnts_ext_mesh, xadj, adjncy, nnodes_elmnts, nodes_elmnts,1)
+
+ print *, 'ZZZZ'
+
+ elmnts_ext_mesh(:,:) = elmnts_ext_mesh(:,:) + 1
+ adjncy(:) = adjncy(:) + 1
+ xadj(:) = xadj(:) + 1
+
+ allocate(elmnts_ext_mesh_sub(ESIZE,nelmnts_ext_mesh*NSUB*NSUB*NSUB))
+ allocate(nodes_coords_ext_mesh_sub(NDIM,ESIZE*nelmnts_ext_mesh*(NSUB+1)*(NSUB+1)*(NSUB+1)))
+ allocate(mat_ext_mesh_sub(nelmnts_ext_mesh*NSUB*NSUB*NSUB))
+
+ nnodes_ext_mesh_sub = 0
+
+ do ispec = 1, nelmnts_ext_mesh
+
+ do ix = 1, NSUB+1
+
+ temporary_nodes(1,ix,1,1) = nodes_coords_ext_mesh(1,elmnts_ext_mesh(1,ispec)) + &
+ ( (nodes_coords_ext_mesh(1,elmnts_ext_mesh(2,ispec)) - nodes_coords_ext_mesh(1,elmnts_ext_mesh(1,ispec))) &
+ / real(NSUB)) * (ix-1)
+ temporary_nodes(2,ix,1,1) = nodes_coords_ext_mesh(2,elmnts_ext_mesh(1,ispec)) + &
+ ( (nodes_coords_ext_mesh(2,elmnts_ext_mesh(2,ispec)) - nodes_coords_ext_mesh(2,elmnts_ext_mesh(1,ispec))) &
+ / real(NSUB)) * (ix-1)
+ temporary_nodes(3,ix,1,1) = nodes_coords_ext_mesh(3,elmnts_ext_mesh(1,ispec)) + &
+ ( (nodes_coords_ext_mesh(3,elmnts_ext_mesh(2,ispec)) - nodes_coords_ext_mesh(3,elmnts_ext_mesh(1,ispec))) &
+ / real(NSUB)) * (ix-1)
+
+ temporary_nodes(1,ix,NSUB+1,1) = nodes_coords_ext_mesh(1,elmnts_ext_mesh(4,ispec)) + &
+ ( (nodes_coords_ext_mesh(1,elmnts_ext_mesh(3,ispec)) - nodes_coords_ext_mesh(1,elmnts_ext_mesh(4,ispec))) &
+ / real(NSUB)) * (ix-1)
+ temporary_nodes(2,ix,NSUB+1,1) = nodes_coords_ext_mesh(2,elmnts_ext_mesh(4,ispec)) + &
+ ( (nodes_coords_ext_mesh(2,elmnts_ext_mesh(3,ispec)) - nodes_coords_ext_mesh(2,elmnts_ext_mesh(4,ispec))) &
+ / real(NSUB)) * (ix-1)
+ temporary_nodes(3,ix,NSUB+1,1) = nodes_coords_ext_mesh(3,elmnts_ext_mesh(4,ispec)) + &
+ ( (nodes_coords_ext_mesh(3,elmnts_ext_mesh(3,ispec)) - nodes_coords_ext_mesh(3,elmnts_ext_mesh(4,ispec))) &
+ / real(NSUB)) * (ix-1)
+
+ temporary_nodes(1,ix,1,NSUB+1) = nodes_coords_ext_mesh(1,elmnts_ext_mesh(5,ispec)) + &
+ ( (nodes_coords_ext_mesh(1,elmnts_ext_mesh(6,ispec)) - nodes_coords_ext_mesh(1,elmnts_ext_mesh(5,ispec))) &
+ / real(NSUB)) * (ix-1)
+ temporary_nodes(2,ix,1,NSUB+1) = nodes_coords_ext_mesh(2,elmnts_ext_mesh(5,ispec)) + &
+ ( (nodes_coords_ext_mesh(2,elmnts_ext_mesh(6,ispec)) - nodes_coords_ext_mesh(2,elmnts_ext_mesh(5,ispec))) &
+ / real(NSUB)) * (ix-1)
+ temporary_nodes(3,ix,1,NSUB+1) = nodes_coords_ext_mesh(3,elmnts_ext_mesh(5,ispec)) + &
+ ( (nodes_coords_ext_mesh(3,elmnts_ext_mesh(6,ispec)) - nodes_coords_ext_mesh(3,elmnts_ext_mesh(5,ispec))) &
+ / real(NSUB)) * (ix-1)
+
+ temporary_nodes(1,ix,NSUB+1,NSUB+1) = nodes_coords_ext_mesh(1,elmnts_ext_mesh(8,ispec)) + &
+ ( (nodes_coords_ext_mesh(1,elmnts_ext_mesh(7,ispec)) - nodes_coords_ext_mesh(1,elmnts_ext_mesh(8,ispec))) &
+ / real(NSUB)) * (ix-1)
+ temporary_nodes(2,ix,NSUB+1,NSUB+1) = nodes_coords_ext_mesh(2,elmnts_ext_mesh(8,ispec)) + &
+ ( (nodes_coords_ext_mesh(2,elmnts_ext_mesh(7,ispec)) - nodes_coords_ext_mesh(2,elmnts_ext_mesh(8,ispec))) &
+ / real(NSUB)) * (ix-1)
+ temporary_nodes(3,ix,NSUB+1,NSUB+1) = nodes_coords_ext_mesh(3,elmnts_ext_mesh(8,ispec)) + &
+ ( (nodes_coords_ext_mesh(3,elmnts_ext_mesh(7,ispec)) - nodes_coords_ext_mesh(3,elmnts_ext_mesh(8,ispec))) &
+ / real(NSUB)) * (ix-1)
+
+ enddo
+
+ do iy = 1, NSUB+1
+
+ temporary_nodes(1,1,iy,1) = nodes_coords_ext_mesh(1,elmnts_ext_mesh(1,ispec)) + &
+ ( (nodes_coords_ext_mesh(1,elmnts_ext_mesh(4,ispec)) - nodes_coords_ext_mesh(1,elmnts_ext_mesh(1,ispec))) &
+ / real(NSUB)) * (iy-1)
+ temporary_nodes(2,1,iy,1) = nodes_coords_ext_mesh(2,elmnts_ext_mesh(1,ispec)) + &
+ ( (nodes_coords_ext_mesh(2,elmnts_ext_mesh(4,ispec)) - nodes_coords_ext_mesh(2,elmnts_ext_mesh(1,ispec))) &
+ / real(NSUB)) * (iy-1)
+ temporary_nodes(3,1,iy,1) = nodes_coords_ext_mesh(3,elmnts_ext_mesh(1,ispec)) + &
+ ( (nodes_coords_ext_mesh(3,elmnts_ext_mesh(4,ispec)) - nodes_coords_ext_mesh(3,elmnts_ext_mesh(1,ispec))) &
+ / real(NSUB)) * (iy-1)
+
+ temporary_nodes(1,NSUB+1,iy,1) = nodes_coords_ext_mesh(1,elmnts_ext_mesh(2,ispec)) + &
+ ( (nodes_coords_ext_mesh(1,elmnts_ext_mesh(3,ispec)) - nodes_coords_ext_mesh(1,elmnts_ext_mesh(2,ispec))) &
+ / real(NSUB)) * (iy-1)
+ temporary_nodes(2,NSUB+1,iy,1) = nodes_coords_ext_mesh(2,elmnts_ext_mesh(2,ispec)) + &
+ ( (nodes_coords_ext_mesh(2,elmnts_ext_mesh(3,ispec)) - nodes_coords_ext_mesh(2,elmnts_ext_mesh(2,ispec))) &
+ / real(NSUB)) * (iy-1)
+ temporary_nodes(3,NSUB+1,iy,1) = nodes_coords_ext_mesh(3,elmnts_ext_mesh(2,ispec)) + &
+ ( (nodes_coords_ext_mesh(3,elmnts_ext_mesh(3,ispec)) - nodes_coords_ext_mesh(3,elmnts_ext_mesh(2,ispec))) &
+ / real(NSUB)) * (iy-1)
+
+ temporary_nodes(1,1,iy,NSUB+1) = nodes_coords_ext_mesh(1,elmnts_ext_mesh(5,ispec)) + &
+ ( (nodes_coords_ext_mesh(1,elmnts_ext_mesh(8,ispec)) - nodes_coords_ext_mesh(1,elmnts_ext_mesh(5,ispec))) &
+ / real(NSUB)) * (iy-1)
+ temporary_nodes(2,1,iy,NSUB+1) = nodes_coords_ext_mesh(2,elmnts_ext_mesh(5,ispec)) + &
+ ( (nodes_coords_ext_mesh(2,elmnts_ext_mesh(8,ispec)) - nodes_coords_ext_mesh(2,elmnts_ext_mesh(5,ispec))) &
+ / real(NSUB)) * (iy-1)
+ temporary_nodes(3,1,iy,NSUB+1) = nodes_coords_ext_mesh(3,elmnts_ext_mesh(5,ispec)) + &
+ ( (nodes_coords_ext_mesh(3,elmnts_ext_mesh(8,ispec)) - nodes_coords_ext_mesh(3,elmnts_ext_mesh(5,ispec))) &
+ / real(NSUB)) * (iy-1)
+
+ temporary_nodes(1,NSUB+1,iy,NSUB+1) = nodes_coords_ext_mesh(1,elmnts_ext_mesh(6,ispec)) + &
+ ( (nodes_coords_ext_mesh(1,elmnts_ext_mesh(7,ispec)) - nodes_coords_ext_mesh(1,elmnts_ext_mesh(6,ispec))) &
+ / real(NSUB)) * (iy-1)
+ temporary_nodes(2,NSUB+1,iy,NSUB+1) = nodes_coords_ext_mesh(2,elmnts_ext_mesh(6,ispec)) + &
+ ( (nodes_coords_ext_mesh(2,elmnts_ext_mesh(7,ispec)) - nodes_coords_ext_mesh(2,elmnts_ext_mesh(6,ispec))) &
+ / real(NSUB)) * (iy-1)
+ temporary_nodes(3,NSUB+1,iy,NSUB+1) = nodes_coords_ext_mesh(3,elmnts_ext_mesh(6,ispec)) + &
+ ( (nodes_coords_ext_mesh(3,elmnts_ext_mesh(7,ispec)) - nodes_coords_ext_mesh(3,elmnts_ext_mesh(6,ispec))) &
+ / real(NSUB)) * (iy-1)
+
+ enddo
+
+ do iz = 1, NSUB+1
+
+ temporary_nodes(1,1,1,iz) = nodes_coords_ext_mesh(1,elmnts_ext_mesh(1,ispec)) + &
+ ( (nodes_coords_ext_mesh(1,elmnts_ext_mesh(5,ispec)) - nodes_coords_ext_mesh(1,elmnts_ext_mesh(1,ispec))) &
+ / real(NSUB)) * (iz-1)
+ temporary_nodes(2,1,1,iz) = nodes_coords_ext_mesh(2,elmnts_ext_mesh(1,ispec)) + &
+ ( (nodes_coords_ext_mesh(2,elmnts_ext_mesh(5,ispec)) - nodes_coords_ext_mesh(2,elmnts_ext_mesh(1,ispec))) &
+ / real(NSUB)) * (iz-1)
+ temporary_nodes(3,1,1,iz) = nodes_coords_ext_mesh(3,elmnts_ext_mesh(1,ispec)) + &
+ ( (nodes_coords_ext_mesh(3,elmnts_ext_mesh(5,ispec)) - nodes_coords_ext_mesh(3,elmnts_ext_mesh(1,ispec))) &
+ / real(NSUB)) * (iz-1)
+
+ temporary_nodes(1,NSUB+1,1,iz) = nodes_coords_ext_mesh(1,elmnts_ext_mesh(2,ispec)) + &
+ ( (nodes_coords_ext_mesh(1,elmnts_ext_mesh(6,ispec)) - nodes_coords_ext_mesh(1,elmnts_ext_mesh(2,ispec))) &
+ / real(NSUB)) * (iz-1)
+ temporary_nodes(2,NSUB+1,1,iz) = nodes_coords_ext_mesh(2,elmnts_ext_mesh(2,ispec)) + &
+ ( (nodes_coords_ext_mesh(2,elmnts_ext_mesh(6,ispec)) - nodes_coords_ext_mesh(2,elmnts_ext_mesh(2,ispec))) &
+ / real(NSUB)) * (iz-1)
+ temporary_nodes(3,NSUB+1,1,iz) = nodes_coords_ext_mesh(3,elmnts_ext_mesh(2,ispec)) + &
+ ( (nodes_coords_ext_mesh(3,elmnts_ext_mesh(6,ispec)) - nodes_coords_ext_mesh(3,elmnts_ext_mesh(2,ispec))) &
+ / real(NSUB)) * (iz-1)
+
+ temporary_nodes(1,1,NSUB+1,iz) = nodes_coords_ext_mesh(1,elmnts_ext_mesh(4,ispec)) + &
+ ( (nodes_coords_ext_mesh(1,elmnts_ext_mesh(8,ispec)) - nodes_coords_ext_mesh(1,elmnts_ext_mesh(4,ispec))) &
+ / real(NSUB)) * (iz-1)
+ temporary_nodes(2,1,NSUB+1,iz) = nodes_coords_ext_mesh(2,elmnts_ext_mesh(4,ispec)) + &
+ ( (nodes_coords_ext_mesh(2,elmnts_ext_mesh(8,ispec)) - nodes_coords_ext_mesh(2,elmnts_ext_mesh(4,ispec))) &
+ / real(NSUB)) * (iz-1)
+ temporary_nodes(3,1,NSUB+1,iz) = nodes_coords_ext_mesh(3,elmnts_ext_mesh(4,ispec)) + &
+ ( (nodes_coords_ext_mesh(3,elmnts_ext_mesh(8,ispec)) - nodes_coords_ext_mesh(3,elmnts_ext_mesh(4,ispec))) &
+ / real(NSUB)) * (iz-1)
+
+ temporary_nodes(1,NSUB+1,NSUB+1,iz) = nodes_coords_ext_mesh(1,elmnts_ext_mesh(3,ispec)) + &
+ ( (nodes_coords_ext_mesh(1,elmnts_ext_mesh(7,ispec)) - nodes_coords_ext_mesh(1,elmnts_ext_mesh(3,ispec))) &
+ / real(NSUB)) * (iz-1)
+ temporary_nodes(2,NSUB+1,NSUB+1,iz) = nodes_coords_ext_mesh(2,elmnts_ext_mesh(3,ispec)) + &
+ ( (nodes_coords_ext_mesh(2,elmnts_ext_mesh(7,ispec)) - nodes_coords_ext_mesh(2,elmnts_ext_mesh(3,ispec))) &
+ / real(NSUB)) * (iz-1)
+ temporary_nodes(3,NSUB+1,NSUB+1,iz) = nodes_coords_ext_mesh(3,elmnts_ext_mesh(3,ispec)) + &
+ ( (nodes_coords_ext_mesh(3,elmnts_ext_mesh(7,ispec)) - nodes_coords_ext_mesh(3,elmnts_ext_mesh(3,ispec))) &
+ / real(NSUB)) * (iz-1)
+
+ enddo
+
+ ix = 1
+ do iy = 2, NSUB
+ do iz = 2, NSUB
+ do idim = 1,NDIM
+ temporary_nodes(idim,ix,iy,iz) = ((temporary_nodes(idim,ix,1,iz) + &
+ ((temporary_nodes(idim,ix,NSUB+1,iz)-temporary_nodes(idim,ix,1,iz)) &
+ / real(NSUB)) * (iy-1)) &
+ + (temporary_nodes(idim,ix,iy,1) + ((temporary_nodes(idim,ix,iy,NSUB+1)-temporary_nodes(idim,ix,iy,1)) &
+ / real(NSUB)) * (iz-1))) &
+ * 1./2.
+
+ enddo
+ enddo
+ enddo
+
+ ix = NSUB+1
+ do iy = 2, NSUB
+ do iz = 2, NSUB
+ do idim = 1,NDIM
+ temporary_nodes(idim,ix,iy,iz) = ((temporary_nodes(idim,ix,1,iz) + &
+ ((temporary_nodes(idim,ix,NSUB+1,iz)-temporary_nodes(idim,ix,1,iz)) &
+ / real(NSUB)) * (iy-1)) &
+ + (temporary_nodes(idim,ix,iy,1) + ((temporary_nodes(idim,ix,iy,NSUB+1)-temporary_nodes(idim,ix,iy,1)) &
+ / real(NSUB)) * (iz-1))) &
+ * 1./2.
+
+ enddo
+ enddo
+ enddo
+
+ iy = 1
+ do ix = 2, NSUB
+ do iz = 2, NSUB
+ do idim = 1,NDIM
+ temporary_nodes(idim,ix,iy,iz) = ((temporary_nodes(idim,1,iy,iz) + &
+ ((temporary_nodes(idim,NSUB+1,iy,iz)-temporary_nodes(idim,1,iy,iz)) &
+ / real(NSUB)) * (ix-1)) &
+ + (temporary_nodes(idim,ix,iy,1) + ((temporary_nodes(idim,ix,iy,NSUB+1)-temporary_nodes(idim,ix,iy,1)) &
+ / real(NSUB)) * (iz-1))) &
+ * 1./2.
+
+ enddo
+ enddo
+ enddo
+
+ iy = NSUB+1
+ do ix = 2, NSUB
+ do iz = 2, NSUB
+ do idim = 1,NDIM
+ temporary_nodes(idim,ix,iy,iz) = ((temporary_nodes(idim,1,iy,iz) + &
+ ((temporary_nodes(idim,NSUB+1,iy,iz)-temporary_nodes(idim,1,iy,iz)) &
+ / real(NSUB)) * (ix-1)) &
+ + (temporary_nodes(idim,ix,iy,1) + ((temporary_nodes(idim,ix,iy,NSUB+1)-temporary_nodes(idim,ix,iy,1)) &
+ / real(NSUB)) * (iz-1))) &
+ * 1./2.
+
+ enddo
+ enddo
+ enddo
+
+ iz = 1
+ do ix = 2, NSUB
+ do iy = 2, NSUB
+ do idim = 1,NDIM
+ temporary_nodes(idim,ix,iy,iz) = ((temporary_nodes(idim,1,iy,iz) + &
+ ((temporary_nodes(idim,NSUB+1,iy,iz)-temporary_nodes(idim,1,iy,iz)) &
+ / real(NSUB)) * (ix-1)) &
+ + (temporary_nodes(idim,ix,1,iz) + ((temporary_nodes(idim,ix,NSUB+1,iz)-temporary_nodes(idim,ix,1,iz)) &
+ / real(NSUB)) * (iy-1))) &
+ * 1./2.
+
+ enddo
+ enddo
+ enddo
+
+ iz = NSUB+1
+ do ix = 2, NSUB
+ do iy = 2, NSUB
+ do idim = 1,NDIM
+ temporary_nodes(idim,ix,iy,iz) = ((temporary_nodes(idim,1,iy,iz) + &
+ ((temporary_nodes(idim,NSUB+1,iy,iz)-temporary_nodes(idim,1,iy,iz)) &
+ / real(NSUB)) * (ix-1)) &
+ + (temporary_nodes(idim,ix,1,iz) + ((temporary_nodes(idim,ix,NSUB+1,iz)-temporary_nodes(idim,ix,1,iz)) &
+ / real(NSUB)) * (iy-1))) &
+ * 1./2.
+
+ enddo
+ enddo
+ enddo
+
+ do ix = 2, NSUB
+ do iy = 2, NSUB
+ do iz = 2, NSUB
+ do idim = 1,NDIM
+ temporary_nodes(idim,ix,iy,iz) = ((temporary_nodes(idim,1,iy,iz) + &
+ ((temporary_nodes(idim,NSUB+1,iy,iz)-temporary_nodes(idim,1,iy,iz)) &
+ / real(NSUB)) * (ix-1)) &
+ + (temporary_nodes(idim,ix,1,iz) + ((temporary_nodes(idim,ix,NSUB+1,iz)-temporary_nodes(idim,ix,1,iz)) &
+ / real(NSUB)) * (iy-1)) &
+ + (temporary_nodes(idim,ix,iy,1) + ((temporary_nodes(idim,ix,iy,NSUB+1)-temporary_nodes(idim,ix,iy,1)) &
+ / real(NSUB)) * (iz-1))) &
+ * 1./3.
+
+ enddo
+ enddo
+ enddo
+ enddo
+
+ temporary_nodes_lookup(:,:,:) = 0
+
+ do ispec_neighbours = xadj(ispec), xadj(ispec+1)-1
+ if ( adjncy(ispec_neighbours) < ispec ) then
+ do ispec_neighbours_sub = (adjncy(ispec_neighbours)-1)*NSUB*NSUB*NSUB + 1, adjncy(ispec_neighbours)*NSUB*NSUB*NSUB
+
+ do ix = 1, NSUB+1
+ do iy = 1, NSUB+1
+ do iz = 1, NSUB+1
+ do inode = 1, ESIZE
+ if ( sqrt( &
+ (temporary_nodes(1,ix,iy,iz)-nodes_coords_ext_mesh_sub(1,elmnts_ext_mesh_sub(inode,ispec_neighbours_sub)))**2 + &
+ (temporary_nodes(2,ix,iy,iz)-nodes_coords_ext_mesh_sub(2,elmnts_ext_mesh_sub(inode,ispec_neighbours_sub)))**2 + &
+ (temporary_nodes(3,ix,iy,iz)-nodes_coords_ext_mesh_sub(3,elmnts_ext_mesh_sub(inode,ispec_neighbours_sub)))**2 ) &
+ < xtol ) then
+ temporary_nodes_lookup(ix,iy,iz) = elmnts_ext_mesh_sub(inode,ispec_neighbours_sub)
+ end if
+
+ enddo
+ enddo
+ enddo
+ enddo
+ enddo
+ end if
+ enddo
+
+ do ix = 1, NSUB+1
+ do iy = 1, NSUB+1
+ do iz = 1, NSUB+1
+ if (temporary_nodes_lookup(ix,iy,iz) == 0 ) then
+ nnodes_ext_mesh_sub = nnodes_ext_mesh_sub + 1
+ temporary_nodes_lookup(ix,iy,iz) = nnodes_ext_mesh_sub
+ nodes_coords_ext_mesh_sub(1,nnodes_ext_mesh_sub) = temporary_nodes(1,ix,iy,iz)
+ nodes_coords_ext_mesh_sub(2,nnodes_ext_mesh_sub) = temporary_nodes(2,ix,iy,iz)
+ nodes_coords_ext_mesh_sub(3,nnodes_ext_mesh_sub) = temporary_nodes(3,ix,iy,iz)
+ end if
+ enddo
+ enddo
+ enddo
+
+ do ix = 1, NSUB
+ do iy = 1, NSUB
+ do iz = 1, NSUB
+ elmnts_ext_mesh_sub(1,(ispec-1)*NSUB*NSUB*NSUB+(ix-1)*NSUB*NSUB+(iy-1)*NSUB+iz) = temporary_nodes_lookup(ix,iy,iz)
+ elmnts_ext_mesh_sub(2,(ispec-1)*NSUB*NSUB*NSUB+(ix-1)*NSUB*NSUB+(iy-1)*NSUB+iz) = temporary_nodes_lookup(ix+1,iy,iz)
+ elmnts_ext_mesh_sub(3,(ispec-1)*NSUB*NSUB*NSUB+(ix-1)*NSUB*NSUB+(iy-1)*NSUB+iz) = temporary_nodes_lookup(ix+1,iy+1,iz)
+ elmnts_ext_mesh_sub(4,(ispec-1)*NSUB*NSUB*NSUB+(ix-1)*NSUB*NSUB+(iy-1)*NSUB+iz) = temporary_nodes_lookup(ix,iy+1,iz)
+ elmnts_ext_mesh_sub(5,(ispec-1)*NSUB*NSUB*NSUB+(ix-1)*NSUB*NSUB+(iy-1)*NSUB+iz) = temporary_nodes_lookup(ix,iy,iz+1)
+ elmnts_ext_mesh_sub(6,(ispec-1)*NSUB*NSUB*NSUB+(ix-1)*NSUB*NSUB+(iy-1)*NSUB+iz) = temporary_nodes_lookup(ix+1,iy,iz+1)
+ elmnts_ext_mesh_sub(7,(ispec-1)*NSUB*NSUB*NSUB+(ix-1)*NSUB*NSUB+(iy-1)*NSUB+iz) = temporary_nodes_lookup(ix+1,iy+1,iz+1)
+ elmnts_ext_mesh_sub(8,(ispec-1)*NSUB*NSUB*NSUB+(ix-1)*NSUB*NSUB+(iy-1)*NSUB+iz) = temporary_nodes_lookup(ix,iy+1,iz+1)
+
+ mat_ext_mesh_sub((ispec-1)*NSUB*NSUB*NSUB+(ix-1)*NSUB*NSUB+(iy-1)*NSUB+iz) = mat_ext_mesh(ispec)
+
+ enddo
+ enddo
+ enddo
+
+ enddo
+
+! check that there really are 8 nodes per element.
+ do ispec = 1, nelmnts_ext_mesh*NSUB*NSUB*NSUB
+ do inode = 1, ESIZE
+ do ix = inode+1, ESIZE
+ if (elmnts_ext_mesh_sub(inode,ispec) == elmnts_ext_mesh_sub(ix,ispec)) then
+ stop 'ERRORERROR'
+ endif
+ enddo
+
+ enddo
+enddo
+
+
+ print *, 'xmin', minval(nodes_coords_ext_mesh_sub(1,:))
+ print *, 'xmax', maxval(nodes_coords_ext_mesh_sub(1,:))
+ print *, 'ymin', minval(nodes_coords_ext_mesh_sub(2,:))
+ print *, 'ymax', maxval(nodes_coords_ext_mesh_sub(2,:))
+ print *, 'zmin', minval(nodes_coords_ext_mesh_sub(3,:))
+ print *, 'zmax', maxval(nodes_coords_ext_mesh_sub(3,:))
+
+
+ open(unit=99, file='./mesh_sub', status='unknown', form='formatted')
+ write(99,*) nelmnts_ext_mesh*NSUB*NSUB*NSUB
+ do ispec = 1, nelmnts_ext_mesh*NSUB*NSUB*NSUB
+ write(99,*) elmnts_ext_mesh_sub(1,ispec), elmnts_ext_mesh_sub(2,ispec), elmnts_ext_mesh_sub(3,ispec), &
+ elmnts_ext_mesh_sub(4,ispec), elmnts_ext_mesh_sub(5,ispec), elmnts_ext_mesh_sub(6,ispec), &
+ elmnts_ext_mesh_sub(7,ispec), elmnts_ext_mesh_sub(8,ispec)
+ end do
+ close(99)
+
+ open(unit=99, file='./mat_sub', status='unknown', form='formatted')
+ do ispec = 1, nelmnts_ext_mesh*NSUB*NSUB*NSUB
+ write(99,*) mat_ext_mesh_sub(ispec)
+ end do
+ close(99)
+
+
+ open(unit=99, file='./nodes_coords_sub', status='unknown', form='formatted')
+ write(99,*) nnodes_ext_mesh_sub
+ do inode = 1, nnodes_ext_mesh_sub
+ write(99,*) nodes_coords_ext_mesh_sub(1,inode), nodes_coords_ext_mesh_sub(2,inode), nodes_coords_ext_mesh_sub(3,inode)
+ end do
+ close(99)
+
+ end program decimate_mesh
+
+
+ !-----------------------------------------------
+ ! Creating dual graph (adjacency is defined by 'ncommonnodes' between two elements).
+ !-----------------------------------------------
+ subroutine mesh2dual_ncommonnodes(nelmnts, nnodes, elmnts, xadj, adjncy, nnodes_elmnts, nodes_elmnts, ncommonnodes)
+
+ include './constants.h'
+
+ integer, intent(in) :: nelmnts
+ integer, intent(in) :: nnodes
+ integer, dimension(0:esize*nelmnts-1), intent(in) :: elmnts
+ integer, dimension(0:nelmnts) :: xadj
+ integer, dimension(0:MAX_NEIGHBOURS*nelmnts-1) :: adjncy
+ integer, dimension(0:nnodes-1) :: nnodes_elmnts
+ integer, dimension(0:nsize*nnodes-1) :: nodes_elmnts
+ integer, intent(in) :: ncommonnodes
+
+ integer :: i, j, k, l, m, nb_edges
+ logical :: is_neighbour
+ integer :: num_node, n
+ integer :: elem_base, elem_target
+ integer :: connectivity
+
+ print *, 'RRRRRRRRRR'
+
+ !allocate(xadj(0:nelmnts))
+ xadj(:) = 0
+ !allocate(adjncy(0:MAX_NEIGHBOURS*nelmnts-1))
+ adjncy(:) = 0
+ !allocate(nnodes_elmnts(0:nnodes-1))
+ nnodes_elmnts(:) = 0
+ !allocate(nodes_elmnts(0:nsize*nnodes-1))
+ 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
+
+ print *, 'nnodes_elmnts'
+
+ ! 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)*MAX_NEIGHBOURS+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)*MAX_NEIGHBOURS+xadj(nodes_elmnts(k+j*nsize))) = nodes_elmnts(l+j*nsize)
+ xadj(nodes_elmnts(k+j*nsize)) = xadj(nodes_elmnts(k+j*nsize)) + 1
+ adjncy(nodes_elmnts(l+j*nsize)*MAX_NEIGHBOURS+xadj(nodes_elmnts(l+j*nsize))) = nodes_elmnts(k+j*nsize)
+ xadj(nodes_elmnts(l+j*nsize)) = xadj(nodes_elmnts(l+j*nsize)) + 1
+ end if
+ end if
+ end do
+ end do
+ end do
+
+ ! 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*MAX_NEIGHBOURS+j)
+ nb_edges = nb_edges + 1
+ end do
+ end do
+
+ xadj(nelmnts) = nb_edges
+
+
+ end subroutine mesh2dual_ncommonnodes
+
Deleted: seismo/3D/SPECFEM3D_SESAME/trunk/UTILS/external_mesh/decimate_mesh/sub.f90
===================================================================
--- seismo/3D/SPECFEM3D_SESAME/trunk/UTILS/external_mesh/model_asteroid_subdivide/sub.f90 2009-03-12 16:31:24 UTC (rev 14300)
+++ seismo/3D/SPECFEM3D_SESAME/trunk/UTILS/external_mesh/decimate_mesh/sub.f90 2009-03-12 17:42:51 UTC (rev 14302)
@@ -1,593 +0,0 @@
-program subdivide_mesh
-
- implicit none
-
- include './constants.h'
-
- integer :: nelmnts_ext_mesh
- integer, dimension(:,:), allocatable :: elmnts_ext_mesh
- integer, dimension(:,:), allocatable :: elmnts_ext_mesh_sub
- integer, dimension(:), allocatable :: mat_ext_mesh
- integer, dimension(:), allocatable :: mat_ext_mesh_sub
-
- integer :: nnodes_ext_mesh
- real, dimension(:,:), allocatable :: nodes_coords_ext_mesh
- real, dimension(:,:), allocatable :: nodes_coords_ext_mesh_sub
-
-
-
- real, dimension(NDIM,NSUB+1,NSUB+1,NSUB+1) :: temporary_nodes
- integer, dimension(NSUB+1,NSUB+1,NSUB+1) :: temporary_nodes_lookup
-
- integer, dimension(:), allocatable :: xadj
- integer, dimension(:), allocatable :: adjncy
- integer, dimension(:), allocatable :: nnodes_elmnts
- integer, dimension(:), allocatable :: nodes_elmnts
-
- integer :: ispec, inode, ispec_neighbours, ispec_neighbours_sub
- integer :: nnodes_ext_mesh_sub
- integer :: i, j, k
- integer :: ix, iy, iz
- integer :: idim
-
- real :: xtol
-
- real :: xminval,yminval,xmaxval,ymaxval,xtypdist,zminval,zmaxval
-
-
-
- open(unit=98, file='./mesh', status='old', form='formatted')
- read(98,*) nelmnts_ext_mesh
- allocate(elmnts_ext_mesh(ESIZE,nelmnts_ext_mesh))
- do ispec = 1, nelmnts_ext_mesh
- read(98,*) elmnts_ext_mesh(1,ispec), elmnts_ext_mesh(2,ispec), elmnts_ext_mesh(3,ispec), elmnts_ext_mesh(4,ispec), &
- elmnts_ext_mesh(5,ispec), elmnts_ext_mesh(6,ispec), elmnts_ext_mesh(7,ispec), elmnts_ext_mesh(8,ispec)
- end do
- close(98)
-
- open(unit=98, file='./mat', status='old', form='formatted')
- allocate(mat_ext_mesh(nelmnts_ext_mesh))
- do ispec = 1, nelmnts_ext_mesh
- read(98,*) mat_ext_mesh(ispec)
- end do
- close(98)
-
- open(unit=98, file='./nodes_coords', status='old', form='formatted')
- read(98,*) nnodes_ext_mesh
- allocate(nodes_coords_ext_mesh(NDIM,nnodes_ext_mesh))
- do inode = 1, nnodes_ext_mesh
- read(98,*) nodes_coords_ext_mesh(1,inode), nodes_coords_ext_mesh(2,inode),nodes_coords_ext_mesh(3,inode)
- end do
- close(98)
-
-
-! check that there really are 8 nodes per element.
- do ispec = 1, nelmnts_ext_mesh
- do inode = 1, ESIZE
- do ix = inode+1, ESIZE
- if (elmnts_ext_mesh(inode,ispec) == elmnts_ext_mesh(ix,ispec)) then
- stop 'ERRORERROR'
- endif
- enddo
-
- enddo
-enddo
-
-! set up local geometric tolerances
- xtypdist=+HUGEVAL
-
- do ispec=1,nelmnts_ext_mesh
-
- xminval=+HUGEVAL
- yminval=+HUGEVAL
- zminval=+HUGEVAL
- xmaxval=-HUGEVAL
- ymaxval=-HUGEVAL
- zmaxval=-HUGEVAL
-
- do inode = 1, 8
- xmaxval=max(nodes_coords_ext_mesh(1,elmnts_ext_mesh(inode,ispec)),xmaxval)
- xminval=min(nodes_coords_ext_mesh(1,elmnts_ext_mesh(inode,ispec)),xminval)
- ymaxval=max(nodes_coords_ext_mesh(2,elmnts_ext_mesh(inode,ispec)),ymaxval)
- yminval=min(nodes_coords_ext_mesh(2,elmnts_ext_mesh(inode,ispec)),yminval)
- zmaxval=max(nodes_coords_ext_mesh(3,elmnts_ext_mesh(inode,ispec)),zmaxval)
- zminval=min(nodes_coords_ext_mesh(3,elmnts_ext_mesh(inode,ispec)),zminval)
- enddo
-
-! compute the minimum typical "size" of an element in the mesh
- xtypdist = min(xtypdist,xmaxval-xminval)
- xtypdist = min(xtypdist,ymaxval-yminval)
- xtypdist = min(xtypdist,zmaxval-zminval)
- !xtypdist = min(xtypdist,sqrt((xmaxval-xminval)**2 + (ymaxval-yminval)**2 + (zmaxval-zminval)**2))
-
- enddo
-
-! define a tolerance, small with respect to the minimum size
- xtol=smallval_tol*xtypdist*1.d7
-
- print *, 'xtypdist' , xtypdist
- print *, 'facteur de tolerance XTOL = ', xtol
-
- print *, 'xmin', minval(nodes_coords_ext_mesh(1,:))
- print *, 'xmax', maxval(nodes_coords_ext_mesh(1,:))
- print *, 'ymin', minval(nodes_coords_ext_mesh(2,:))
- print *, 'ymax', maxval(nodes_coords_ext_mesh(2,:))
- print *, 'zmin', minval(nodes_coords_ext_mesh(3,:))
- print *, 'zmax', maxval(nodes_coords_ext_mesh(3,:))
-
-
-
-! we build the graph
- elmnts_ext_mesh(:,:) = elmnts_ext_mesh(:,:) - 1
-
- allocate(xadj(1:nelmnts_ext_mesh+1))
- allocate(adjncy(1:MAX_NEIGHBOURS*nelmnts_ext_mesh))
- allocate(nnodes_elmnts(1:nnodes_ext_mesh))
- allocate(nodes_elmnts(1:NSIZE*nnodes_ext_mesh))
-
- call mesh2dual_ncommonnodes(nelmnts_ext_mesh, nnodes_ext_mesh, elmnts_ext_mesh, xadj, adjncy, nnodes_elmnts, nodes_elmnts,1)
-
- print *, 'ZZZZ'
-
- elmnts_ext_mesh(:,:) = elmnts_ext_mesh(:,:) + 1
- adjncy(:) = adjncy(:) + 1
- xadj(:) = xadj(:) + 1
-
- allocate(elmnts_ext_mesh_sub(ESIZE,nelmnts_ext_mesh*NSUB*NSUB*NSUB))
- allocate(nodes_coords_ext_mesh_sub(NDIM,ESIZE*nelmnts_ext_mesh*(NSUB+1)*(NSUB+1)*(NSUB+1)))
- allocate(mat_ext_mesh_sub(nelmnts_ext_mesh*NSUB*NSUB*NSUB))
-
- nnodes_ext_mesh_sub = 0
-
- do ispec = 1, nelmnts_ext_mesh
-
- do ix = 1, NSUB+1
-
- temporary_nodes(1,ix,1,1) = nodes_coords_ext_mesh(1,elmnts_ext_mesh(1,ispec)) + &
- ( (nodes_coords_ext_mesh(1,elmnts_ext_mesh(2,ispec)) - nodes_coords_ext_mesh(1,elmnts_ext_mesh(1,ispec))) &
- / real(NSUB)) * (ix-1)
- temporary_nodes(2,ix,1,1) = nodes_coords_ext_mesh(2,elmnts_ext_mesh(1,ispec)) + &
- ( (nodes_coords_ext_mesh(2,elmnts_ext_mesh(2,ispec)) - nodes_coords_ext_mesh(2,elmnts_ext_mesh(1,ispec))) &
- / real(NSUB)) * (ix-1)
- temporary_nodes(3,ix,1,1) = nodes_coords_ext_mesh(3,elmnts_ext_mesh(1,ispec)) + &
- ( (nodes_coords_ext_mesh(3,elmnts_ext_mesh(2,ispec)) - nodes_coords_ext_mesh(3,elmnts_ext_mesh(1,ispec))) &
- / real(NSUB)) * (ix-1)
-
- temporary_nodes(1,ix,NSUB+1,1) = nodes_coords_ext_mesh(1,elmnts_ext_mesh(4,ispec)) + &
- ( (nodes_coords_ext_mesh(1,elmnts_ext_mesh(3,ispec)) - nodes_coords_ext_mesh(1,elmnts_ext_mesh(4,ispec))) &
- / real(NSUB)) * (ix-1)
- temporary_nodes(2,ix,NSUB+1,1) = nodes_coords_ext_mesh(2,elmnts_ext_mesh(4,ispec)) + &
- ( (nodes_coords_ext_mesh(2,elmnts_ext_mesh(3,ispec)) - nodes_coords_ext_mesh(2,elmnts_ext_mesh(4,ispec))) &
- / real(NSUB)) * (ix-1)
- temporary_nodes(3,ix,NSUB+1,1) = nodes_coords_ext_mesh(3,elmnts_ext_mesh(4,ispec)) + &
- ( (nodes_coords_ext_mesh(3,elmnts_ext_mesh(3,ispec)) - nodes_coords_ext_mesh(3,elmnts_ext_mesh(4,ispec))) &
- / real(NSUB)) * (ix-1)
-
- temporary_nodes(1,ix,1,NSUB+1) = nodes_coords_ext_mesh(1,elmnts_ext_mesh(5,ispec)) + &
- ( (nodes_coords_ext_mesh(1,elmnts_ext_mesh(6,ispec)) - nodes_coords_ext_mesh(1,elmnts_ext_mesh(5,ispec))) &
- / real(NSUB)) * (ix-1)
- temporary_nodes(2,ix,1,NSUB+1) = nodes_coords_ext_mesh(2,elmnts_ext_mesh(5,ispec)) + &
- ( (nodes_coords_ext_mesh(2,elmnts_ext_mesh(6,ispec)) - nodes_coords_ext_mesh(2,elmnts_ext_mesh(5,ispec))) &
- / real(NSUB)) * (ix-1)
- temporary_nodes(3,ix,1,NSUB+1) = nodes_coords_ext_mesh(3,elmnts_ext_mesh(5,ispec)) + &
- ( (nodes_coords_ext_mesh(3,elmnts_ext_mesh(6,ispec)) - nodes_coords_ext_mesh(3,elmnts_ext_mesh(5,ispec))) &
- / real(NSUB)) * (ix-1)
-
- temporary_nodes(1,ix,NSUB+1,NSUB+1) = nodes_coords_ext_mesh(1,elmnts_ext_mesh(8,ispec)) + &
- ( (nodes_coords_ext_mesh(1,elmnts_ext_mesh(7,ispec)) - nodes_coords_ext_mesh(1,elmnts_ext_mesh(8,ispec))) &
- / real(NSUB)) * (ix-1)
- temporary_nodes(2,ix,NSUB+1,NSUB+1) = nodes_coords_ext_mesh(2,elmnts_ext_mesh(8,ispec)) + &
- ( (nodes_coords_ext_mesh(2,elmnts_ext_mesh(7,ispec)) - nodes_coords_ext_mesh(2,elmnts_ext_mesh(8,ispec))) &
- / real(NSUB)) * (ix-1)
- temporary_nodes(3,ix,NSUB+1,NSUB+1) = nodes_coords_ext_mesh(3,elmnts_ext_mesh(8,ispec)) + &
- ( (nodes_coords_ext_mesh(3,elmnts_ext_mesh(7,ispec)) - nodes_coords_ext_mesh(3,elmnts_ext_mesh(8,ispec))) &
- / real(NSUB)) * (ix-1)
-
- enddo
-
- do iy = 1, NSUB+1
-
- temporary_nodes(1,1,iy,1) = nodes_coords_ext_mesh(1,elmnts_ext_mesh(1,ispec)) + &
- ( (nodes_coords_ext_mesh(1,elmnts_ext_mesh(4,ispec)) - nodes_coords_ext_mesh(1,elmnts_ext_mesh(1,ispec))) &
- / real(NSUB)) * (iy-1)
- temporary_nodes(2,1,iy,1) = nodes_coords_ext_mesh(2,elmnts_ext_mesh(1,ispec)) + &
- ( (nodes_coords_ext_mesh(2,elmnts_ext_mesh(4,ispec)) - nodes_coords_ext_mesh(2,elmnts_ext_mesh(1,ispec))) &
- / real(NSUB)) * (iy-1)
- temporary_nodes(3,1,iy,1) = nodes_coords_ext_mesh(3,elmnts_ext_mesh(1,ispec)) + &
- ( (nodes_coords_ext_mesh(3,elmnts_ext_mesh(4,ispec)) - nodes_coords_ext_mesh(3,elmnts_ext_mesh(1,ispec))) &
- / real(NSUB)) * (iy-1)
-
- temporary_nodes(1,NSUB+1,iy,1) = nodes_coords_ext_mesh(1,elmnts_ext_mesh(2,ispec)) + &
- ( (nodes_coords_ext_mesh(1,elmnts_ext_mesh(3,ispec)) - nodes_coords_ext_mesh(1,elmnts_ext_mesh(2,ispec))) &
- / real(NSUB)) * (iy-1)
- temporary_nodes(2,NSUB+1,iy,1) = nodes_coords_ext_mesh(2,elmnts_ext_mesh(2,ispec)) + &
- ( (nodes_coords_ext_mesh(2,elmnts_ext_mesh(3,ispec)) - nodes_coords_ext_mesh(2,elmnts_ext_mesh(2,ispec))) &
- / real(NSUB)) * (iy-1)
- temporary_nodes(3,NSUB+1,iy,1) = nodes_coords_ext_mesh(3,elmnts_ext_mesh(2,ispec)) + &
- ( (nodes_coords_ext_mesh(3,elmnts_ext_mesh(3,ispec)) - nodes_coords_ext_mesh(3,elmnts_ext_mesh(2,ispec))) &
- / real(NSUB)) * (iy-1)
-
- temporary_nodes(1,1,iy,NSUB+1) = nodes_coords_ext_mesh(1,elmnts_ext_mesh(5,ispec)) + &
- ( (nodes_coords_ext_mesh(1,elmnts_ext_mesh(8,ispec)) - nodes_coords_ext_mesh(1,elmnts_ext_mesh(5,ispec))) &
- / real(NSUB)) * (iy-1)
- temporary_nodes(2,1,iy,NSUB+1) = nodes_coords_ext_mesh(2,elmnts_ext_mesh(5,ispec)) + &
- ( (nodes_coords_ext_mesh(2,elmnts_ext_mesh(8,ispec)) - nodes_coords_ext_mesh(2,elmnts_ext_mesh(5,ispec))) &
- / real(NSUB)) * (iy-1)
- temporary_nodes(3,1,iy,NSUB+1) = nodes_coords_ext_mesh(3,elmnts_ext_mesh(5,ispec)) + &
- ( (nodes_coords_ext_mesh(3,elmnts_ext_mesh(8,ispec)) - nodes_coords_ext_mesh(3,elmnts_ext_mesh(5,ispec))) &
- / real(NSUB)) * (iy-1)
-
- temporary_nodes(1,NSUB+1,iy,NSUB+1) = nodes_coords_ext_mesh(1,elmnts_ext_mesh(6,ispec)) + &
- ( (nodes_coords_ext_mesh(1,elmnts_ext_mesh(7,ispec)) - nodes_coords_ext_mesh(1,elmnts_ext_mesh(6,ispec))) &
- / real(NSUB)) * (iy-1)
- temporary_nodes(2,NSUB+1,iy,NSUB+1) = nodes_coords_ext_mesh(2,elmnts_ext_mesh(6,ispec)) + &
- ( (nodes_coords_ext_mesh(2,elmnts_ext_mesh(7,ispec)) - nodes_coords_ext_mesh(2,elmnts_ext_mesh(6,ispec))) &
- / real(NSUB)) * (iy-1)
- temporary_nodes(3,NSUB+1,iy,NSUB+1) = nodes_coords_ext_mesh(3,elmnts_ext_mesh(6,ispec)) + &
- ( (nodes_coords_ext_mesh(3,elmnts_ext_mesh(7,ispec)) - nodes_coords_ext_mesh(3,elmnts_ext_mesh(6,ispec))) &
- / real(NSUB)) * (iy-1)
-
- enddo
-
- do iz = 1, NSUB+1
-
- temporary_nodes(1,1,1,iz) = nodes_coords_ext_mesh(1,elmnts_ext_mesh(1,ispec)) + &
- ( (nodes_coords_ext_mesh(1,elmnts_ext_mesh(5,ispec)) - nodes_coords_ext_mesh(1,elmnts_ext_mesh(1,ispec))) &
- / real(NSUB)) * (iz-1)
- temporary_nodes(2,1,1,iz) = nodes_coords_ext_mesh(2,elmnts_ext_mesh(1,ispec)) + &
- ( (nodes_coords_ext_mesh(2,elmnts_ext_mesh(5,ispec)) - nodes_coords_ext_mesh(2,elmnts_ext_mesh(1,ispec))) &
- / real(NSUB)) * (iz-1)
- temporary_nodes(3,1,1,iz) = nodes_coords_ext_mesh(3,elmnts_ext_mesh(1,ispec)) + &
- ( (nodes_coords_ext_mesh(3,elmnts_ext_mesh(5,ispec)) - nodes_coords_ext_mesh(3,elmnts_ext_mesh(1,ispec))) &
- / real(NSUB)) * (iz-1)
-
- temporary_nodes(1,NSUB+1,1,iz) = nodes_coords_ext_mesh(1,elmnts_ext_mesh(2,ispec)) + &
- ( (nodes_coords_ext_mesh(1,elmnts_ext_mesh(6,ispec)) - nodes_coords_ext_mesh(1,elmnts_ext_mesh(2,ispec))) &
- / real(NSUB)) * (iz-1)
- temporary_nodes(2,NSUB+1,1,iz) = nodes_coords_ext_mesh(2,elmnts_ext_mesh(2,ispec)) + &
- ( (nodes_coords_ext_mesh(2,elmnts_ext_mesh(6,ispec)) - nodes_coords_ext_mesh(2,elmnts_ext_mesh(2,ispec))) &
- / real(NSUB)) * (iz-1)
- temporary_nodes(3,NSUB+1,1,iz) = nodes_coords_ext_mesh(3,elmnts_ext_mesh(2,ispec)) + &
- ( (nodes_coords_ext_mesh(3,elmnts_ext_mesh(6,ispec)) - nodes_coords_ext_mesh(3,elmnts_ext_mesh(2,ispec))) &
- / real(NSUB)) * (iz-1)
-
- temporary_nodes(1,1,NSUB+1,iz) = nodes_coords_ext_mesh(1,elmnts_ext_mesh(4,ispec)) + &
- ( (nodes_coords_ext_mesh(1,elmnts_ext_mesh(8,ispec)) - nodes_coords_ext_mesh(1,elmnts_ext_mesh(4,ispec))) &
- / real(NSUB)) * (iz-1)
- temporary_nodes(2,1,NSUB+1,iz) = nodes_coords_ext_mesh(2,elmnts_ext_mesh(4,ispec)) + &
- ( (nodes_coords_ext_mesh(2,elmnts_ext_mesh(8,ispec)) - nodes_coords_ext_mesh(2,elmnts_ext_mesh(4,ispec))) &
- / real(NSUB)) * (iz-1)
- temporary_nodes(3,1,NSUB+1,iz) = nodes_coords_ext_mesh(3,elmnts_ext_mesh(4,ispec)) + &
- ( (nodes_coords_ext_mesh(3,elmnts_ext_mesh(8,ispec)) - nodes_coords_ext_mesh(3,elmnts_ext_mesh(4,ispec))) &
- / real(NSUB)) * (iz-1)
-
- temporary_nodes(1,NSUB+1,NSUB+1,iz) = nodes_coords_ext_mesh(1,elmnts_ext_mesh(3,ispec)) + &
- ( (nodes_coords_ext_mesh(1,elmnts_ext_mesh(7,ispec)) - nodes_coords_ext_mesh(1,elmnts_ext_mesh(3,ispec))) &
- / real(NSUB)) * (iz-1)
- temporary_nodes(2,NSUB+1,NSUB+1,iz) = nodes_coords_ext_mesh(2,elmnts_ext_mesh(3,ispec)) + &
- ( (nodes_coords_ext_mesh(2,elmnts_ext_mesh(7,ispec)) - nodes_coords_ext_mesh(2,elmnts_ext_mesh(3,ispec))) &
- / real(NSUB)) * (iz-1)
- temporary_nodes(3,NSUB+1,NSUB+1,iz) = nodes_coords_ext_mesh(3,elmnts_ext_mesh(3,ispec)) + &
- ( (nodes_coords_ext_mesh(3,elmnts_ext_mesh(7,ispec)) - nodes_coords_ext_mesh(3,elmnts_ext_mesh(3,ispec))) &
- / real(NSUB)) * (iz-1)
-
- enddo
-
- ix = 1
- do iy = 2, NSUB
- do iz = 2, NSUB
- do idim = 1,NDIM
- temporary_nodes(idim,ix,iy,iz) = ((temporary_nodes(idim,ix,1,iz) + &
- ((temporary_nodes(idim,ix,NSUB+1,iz)-temporary_nodes(idim,ix,1,iz)) &
- / real(NSUB)) * (iy-1)) &
- + (temporary_nodes(idim,ix,iy,1) + ((temporary_nodes(idim,ix,iy,NSUB+1)-temporary_nodes(idim,ix,iy,1)) &
- / real(NSUB)) * (iz-1))) &
- * 1./2.
-
- enddo
- enddo
- enddo
-
- ix = NSUB+1
- do iy = 2, NSUB
- do iz = 2, NSUB
- do idim = 1,NDIM
- temporary_nodes(idim,ix,iy,iz) = ((temporary_nodes(idim,ix,1,iz) + &
- ((temporary_nodes(idim,ix,NSUB+1,iz)-temporary_nodes(idim,ix,1,iz)) &
- / real(NSUB)) * (iy-1)) &
- + (temporary_nodes(idim,ix,iy,1) + ((temporary_nodes(idim,ix,iy,NSUB+1)-temporary_nodes(idim,ix,iy,1)) &
- / real(NSUB)) * (iz-1))) &
- * 1./2.
-
- enddo
- enddo
- enddo
-
- iy = 1
- do ix = 2, NSUB
- do iz = 2, NSUB
- do idim = 1,NDIM
- temporary_nodes(idim,ix,iy,iz) = ((temporary_nodes(idim,1,iy,iz) + &
- ((temporary_nodes(idim,NSUB+1,iy,iz)-temporary_nodes(idim,1,iy,iz)) &
- / real(NSUB)) * (ix-1)) &
- + (temporary_nodes(idim,ix,iy,1) + ((temporary_nodes(idim,ix,iy,NSUB+1)-temporary_nodes(idim,ix,iy,1)) &
- / real(NSUB)) * (iz-1))) &
- * 1./2.
-
- enddo
- enddo
- enddo
-
- iy = NSUB+1
- do ix = 2, NSUB
- do iz = 2, NSUB
- do idim = 1,NDIM
- temporary_nodes(idim,ix,iy,iz) = ((temporary_nodes(idim,1,iy,iz) + &
- ((temporary_nodes(idim,NSUB+1,iy,iz)-temporary_nodes(idim,1,iy,iz)) &
- / real(NSUB)) * (ix-1)) &
- + (temporary_nodes(idim,ix,iy,1) + ((temporary_nodes(idim,ix,iy,NSUB+1)-temporary_nodes(idim,ix,iy,1)) &
- / real(NSUB)) * (iz-1))) &
- * 1./2.
-
- enddo
- enddo
- enddo
-
- iz = 1
- do ix = 2, NSUB
- do iy = 2, NSUB
- do idim = 1,NDIM
- temporary_nodes(idim,ix,iy,iz) = ((temporary_nodes(idim,1,iy,iz) + &
- ((temporary_nodes(idim,NSUB+1,iy,iz)-temporary_nodes(idim,1,iy,iz)) &
- / real(NSUB)) * (ix-1)) &
- + (temporary_nodes(idim,ix,1,iz) + ((temporary_nodes(idim,ix,NSUB+1,iz)-temporary_nodes(idim,ix,1,iz)) &
- / real(NSUB)) * (iy-1))) &
- * 1./2.
-
- enddo
- enddo
- enddo
-
- iz = NSUB+1
- do ix = 2, NSUB
- do iy = 2, NSUB
- do idim = 1,NDIM
- temporary_nodes(idim,ix,iy,iz) = ((temporary_nodes(idim,1,iy,iz) + &
- ((temporary_nodes(idim,NSUB+1,iy,iz)-temporary_nodes(idim,1,iy,iz)) &
- / real(NSUB)) * (ix-1)) &
- + (temporary_nodes(idim,ix,1,iz) + ((temporary_nodes(idim,ix,NSUB+1,iz)-temporary_nodes(idim,ix,1,iz)) &
- / real(NSUB)) * (iy-1))) &
- * 1./2.
-
- enddo
- enddo
- enddo
-
- do ix = 2, NSUB
- do iy = 2, NSUB
- do iz = 2, NSUB
- do idim = 1,NDIM
- temporary_nodes(idim,ix,iy,iz) = ((temporary_nodes(idim,1,iy,iz) + &
- ((temporary_nodes(idim,NSUB+1,iy,iz)-temporary_nodes(idim,1,iy,iz)) &
- / real(NSUB)) * (ix-1)) &
- + (temporary_nodes(idim,ix,1,iz) + ((temporary_nodes(idim,ix,NSUB+1,iz)-temporary_nodes(idim,ix,1,iz)) &
- / real(NSUB)) * (iy-1)) &
- + (temporary_nodes(idim,ix,iy,1) + ((temporary_nodes(idim,ix,iy,NSUB+1)-temporary_nodes(idim,ix,iy,1)) &
- / real(NSUB)) * (iz-1))) &
- * 1./3.
-
- enddo
- enddo
- enddo
- enddo
-
- temporary_nodes_lookup(:,:,:) = 0
-
- do ispec_neighbours = xadj(ispec), xadj(ispec+1)-1
- if ( adjncy(ispec_neighbours) < ispec ) then
- do ispec_neighbours_sub = (adjncy(ispec_neighbours)-1)*NSUB*NSUB*NSUB + 1, adjncy(ispec_neighbours)*NSUB*NSUB*NSUB
-
- do ix = 1, NSUB+1
- do iy = 1, NSUB+1
- do iz = 1, NSUB+1
- do inode = 1, ESIZE
- if ( sqrt( &
- (temporary_nodes(1,ix,iy,iz)-nodes_coords_ext_mesh_sub(1,elmnts_ext_mesh_sub(inode,ispec_neighbours_sub)))**2 + &
- (temporary_nodes(2,ix,iy,iz)-nodes_coords_ext_mesh_sub(2,elmnts_ext_mesh_sub(inode,ispec_neighbours_sub)))**2 + &
- (temporary_nodes(3,ix,iy,iz)-nodes_coords_ext_mesh_sub(3,elmnts_ext_mesh_sub(inode,ispec_neighbours_sub)))**2 ) &
- < xtol ) then
- temporary_nodes_lookup(ix,iy,iz) = elmnts_ext_mesh_sub(inode,ispec_neighbours_sub)
- end if
-
- enddo
- enddo
- enddo
- enddo
- enddo
- end if
- enddo
-
- do ix = 1, NSUB+1
- do iy = 1, NSUB+1
- do iz = 1, NSUB+1
- if (temporary_nodes_lookup(ix,iy,iz) == 0 ) then
- nnodes_ext_mesh_sub = nnodes_ext_mesh_sub + 1
- temporary_nodes_lookup(ix,iy,iz) = nnodes_ext_mesh_sub
- nodes_coords_ext_mesh_sub(1,nnodes_ext_mesh_sub) = temporary_nodes(1,ix,iy,iz)
- nodes_coords_ext_mesh_sub(2,nnodes_ext_mesh_sub) = temporary_nodes(2,ix,iy,iz)
- nodes_coords_ext_mesh_sub(3,nnodes_ext_mesh_sub) = temporary_nodes(3,ix,iy,iz)
- end if
- enddo
- enddo
- enddo
-
- do ix = 1, NSUB
- do iy = 1, NSUB
- do iz = 1, NSUB
- elmnts_ext_mesh_sub(1,(ispec-1)*NSUB*NSUB*NSUB+(ix-1)*NSUB*NSUB+(iy-1)*NSUB+iz) = temporary_nodes_lookup(ix,iy,iz)
- elmnts_ext_mesh_sub(2,(ispec-1)*NSUB*NSUB*NSUB+(ix-1)*NSUB*NSUB+(iy-1)*NSUB+iz) = temporary_nodes_lookup(ix+1,iy,iz)
- elmnts_ext_mesh_sub(3,(ispec-1)*NSUB*NSUB*NSUB+(ix-1)*NSUB*NSUB+(iy-1)*NSUB+iz) = temporary_nodes_lookup(ix+1,iy+1,iz)
- elmnts_ext_mesh_sub(4,(ispec-1)*NSUB*NSUB*NSUB+(ix-1)*NSUB*NSUB+(iy-1)*NSUB+iz) = temporary_nodes_lookup(ix,iy+1,iz)
- elmnts_ext_mesh_sub(5,(ispec-1)*NSUB*NSUB*NSUB+(ix-1)*NSUB*NSUB+(iy-1)*NSUB+iz) = temporary_nodes_lookup(ix,iy,iz+1)
- elmnts_ext_mesh_sub(6,(ispec-1)*NSUB*NSUB*NSUB+(ix-1)*NSUB*NSUB+(iy-1)*NSUB+iz) = temporary_nodes_lookup(ix+1,iy,iz+1)
- elmnts_ext_mesh_sub(7,(ispec-1)*NSUB*NSUB*NSUB+(ix-1)*NSUB*NSUB+(iy-1)*NSUB+iz) = temporary_nodes_lookup(ix+1,iy+1,iz+1)
- elmnts_ext_mesh_sub(8,(ispec-1)*NSUB*NSUB*NSUB+(ix-1)*NSUB*NSUB+(iy-1)*NSUB+iz) = temporary_nodes_lookup(ix,iy+1,iz+1)
-
- mat_ext_mesh_sub((ispec-1)*NSUB*NSUB*NSUB+(ix-1)*NSUB*NSUB+(iy-1)*NSUB+iz) = mat_ext_mesh(ispec)
-
- enddo
- enddo
- enddo
-
- enddo
-
-! check that there really are 8 nodes per element.
- do ispec = 1, nelmnts_ext_mesh*NSUB*NSUB*NSUB
- do inode = 1, ESIZE
- do ix = inode+1, ESIZE
- if (elmnts_ext_mesh_sub(inode,ispec) == elmnts_ext_mesh_sub(ix,ispec)) then
- stop 'ERRORERROR'
- endif
- enddo
-
- enddo
-enddo
-
-
- print *, 'xmin', minval(nodes_coords_ext_mesh_sub(1,:))
- print *, 'xmax', maxval(nodes_coords_ext_mesh_sub(1,:))
- print *, 'ymin', minval(nodes_coords_ext_mesh_sub(2,:))
- print *, 'ymax', maxval(nodes_coords_ext_mesh_sub(2,:))
- print *, 'zmin', minval(nodes_coords_ext_mesh_sub(3,:))
- print *, 'zmax', maxval(nodes_coords_ext_mesh_sub(3,:))
-
-
- open(unit=99, file='./mesh_sub', status='unknown', form='formatted')
- write(99,*) nelmnts_ext_mesh*NSUB*NSUB*NSUB
- do ispec = 1, nelmnts_ext_mesh*NSUB*NSUB*NSUB
- write(99,*) elmnts_ext_mesh_sub(1,ispec), elmnts_ext_mesh_sub(2,ispec), elmnts_ext_mesh_sub(3,ispec), &
- elmnts_ext_mesh_sub(4,ispec), elmnts_ext_mesh_sub(5,ispec), elmnts_ext_mesh_sub(6,ispec), &
- elmnts_ext_mesh_sub(7,ispec), elmnts_ext_mesh_sub(8,ispec)
- end do
- close(99)
-
- open(unit=99, file='./mat_sub', status='unknown', form='formatted')
- do ispec = 1, nelmnts_ext_mesh*NSUB*NSUB*NSUB
- write(99,*) mat_ext_mesh_sub(ispec)
- end do
- close(99)
-
-
- open(unit=99, file='./nodes_coords_sub', status='unknown', form='formatted')
- write(99,*) nnodes_ext_mesh_sub
- do inode = 1, nnodes_ext_mesh_sub
- write(99,*) nodes_coords_ext_mesh_sub(1,inode), nodes_coords_ext_mesh_sub(2,inode), nodes_coords_ext_mesh_sub(3,inode)
- end do
- close(99)
-
-end program subdivide_mesh
-
-
-
- !-----------------------------------------------
- ! Creating dual graph (adjacency is defined by 'ncommonnodes' between two elements).
- !-----------------------------------------------
- subroutine mesh2dual_ncommonnodes(nelmnts, nnodes, elmnts, xadj, adjncy, nnodes_elmnts, nodes_elmnts, ncommonnodes)
-
- include './constants.h'
-
- integer, intent(in) :: nelmnts
- integer, intent(in) :: nnodes
- integer, dimension(0:esize*nelmnts-1), intent(in) :: elmnts
- integer, dimension(0:nelmnts) :: xadj
- integer, dimension(0:MAX_NEIGHBOURS*nelmnts-1) :: adjncy
- integer, dimension(0:nnodes-1) :: nnodes_elmnts
- integer, dimension(0:nsize*nnodes-1) :: nodes_elmnts
- integer, intent(in) :: ncommonnodes
-
- integer :: i, j, k, l, m, nb_edges
- logical :: is_neighbour
- integer :: num_node, n
- integer :: elem_base, elem_target
- integer :: connectivity
-
- print *, 'RRRRRRRRRR'
-
- !allocate(xadj(0:nelmnts))
- xadj(:) = 0
- !allocate(adjncy(0:MAX_NEIGHBOURS*nelmnts-1))
- adjncy(:) = 0
- !allocate(nnodes_elmnts(0:nnodes-1))
- nnodes_elmnts(:) = 0
- !allocate(nodes_elmnts(0:nsize*nnodes-1))
- 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
-
- print *, 'nnodes_elmnts'
-
- ! 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)*MAX_NEIGHBOURS+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)*MAX_NEIGHBOURS+xadj(nodes_elmnts(k+j*nsize))) = nodes_elmnts(l+j*nsize)
- xadj(nodes_elmnts(k+j*nsize)) = xadj(nodes_elmnts(k+j*nsize)) + 1
- adjncy(nodes_elmnts(l+j*nsize)*MAX_NEIGHBOURS+xadj(nodes_elmnts(l+j*nsize))) = nodes_elmnts(k+j*nsize)
- xadj(nodes_elmnts(l+j*nsize)) = xadj(nodes_elmnts(l+j*nsize)) + 1
- end if
- end if
- end do
- end do
- end do
-
- ! 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*MAX_NEIGHBOURS+j)
- nb_edges = nb_edges + 1
- end do
- end do
-
- xadj(nelmnts) = nb_edges
-
-
- end subroutine mesh2dual_ncommonnodes
-
Copied: seismo/3D/SPECFEM3D_SESAME/trunk/UTILS/external_mesh/decompose_mesh (from rev 14300, seismo/3D/SPECFEM3D_SESAME/trunk/UTILS/external_mesh/pre_meshfem3D_asteroid_subdivide)
Property changes on: seismo/3D/SPECFEM3D_SESAME/trunk/UTILS/external_mesh/decompose_mesh
___________________________________________________________________
Name: svn:mergeinfo
+
Deleted: seismo/3D/SPECFEM3D_SESAME/trunk/UTILS/external_mesh/decompose_mesh/compil_all.sh
===================================================================
--- seismo/3D/SPECFEM3D_SESAME/trunk/UTILS/external_mesh/pre_meshfem3D_asteroid_subdivide/compil_all.sh 2009-03-12 16:31:24 UTC (rev 14300)
+++ seismo/3D/SPECFEM3D_SESAME/trunk/UTILS/external_mesh/decompose_mesh/compil_all.sh 2009-03-12 17:42:51 UTC (rev 14302)
@@ -1,13 +0,0 @@
-#!/bin/sh
-
-#. /opt/intel/fce/10.0.026/bin/ifortvars.sh
-# export FCFLAGS="-g -traceback -implicitnone -warn stderrors -warn truncated_source -warn argument_checking -warn unused -warn declarations -std95 -assume byterecl -C"
-
-rm *.o *.mod ./a.out
-
-gfortran -c part_pre_meshfem3D.f90
-gfortran -c pre_meshfem3D.f90
-#gfortran pre_meshfem3D.o part_pre_meshfem3D.o ~/utils/metis-4.0/libmetis.a
-#gfortran pre_meshfem3D.o part_pre_meshfem3D.o ~/utils/scotch_5.1/lib/libscotchmetis.a ~/utils/scotch_5.1/lib/libscotch.a ~/utils/scotch_5.1/lib/libscotcherr.a
-gfortran pre_meshfem3D.o part_pre_meshfem3D.o ~/utils/scotch_5.1/lib/libscotch.a ~/utils/scotch_5.1/lib/libscotcherr.a
-
Copied: seismo/3D/SPECFEM3D_SESAME/trunk/UTILS/external_mesh/decompose_mesh/compile_all.csh (from rev 14301, seismo/3D/SPECFEM3D_SESAME/trunk/UTILS/external_mesh/pre_meshfem3D_asteroid_subdivide/compile_all.csh)
===================================================================
--- seismo/3D/SPECFEM3D_SESAME/trunk/UTILS/external_mesh/decompose_mesh/compile_all.csh (rev 0)
+++ seismo/3D/SPECFEM3D_SESAME/trunk/UTILS/external_mesh/decompose_mesh/compile_all.csh 2009-03-12 17:42:51 UTC (rev 14302)
@@ -0,0 +1,13 @@
+#!/bin/sh
+
+#. /opt/intel/fce/10.0.026/bin/ifortvars.sh
+# export FCFLAGS="-g -traceback -implicitnone -warn stderrors -warn truncated_source -warn argument_checking -warn unused -warn declarations -std95 -assume byterecl -C"
+
+rm *.o *.mod ./a.out
+
+gfortran -c part_decompose_mesh.f90
+gfortran -c decompose_mesh.f90
+#gfortran decompose_mesh.o part_decompose_mesh.o ~/utils/metis-4.0/libmetis.a
+#gfortran decompose_mesh.o part_decompose_mesh.o ~/utils/scotch_5.1/lib/libscotchmetis.a ~/utils/scotch_5.1/lib/libscotch.a ~/utils/scotch_5.1/lib/libscotcherr.a
+gfortran decompose_mesh.o part_decompose_mesh.o ~/utils/scotch_5.1/lib/libscotch.a ~/utils/scotch_5.1/lib/libscotcherr.a
+
Copied: seismo/3D/SPECFEM3D_SESAME/trunk/UTILS/external_mesh/decompose_mesh/constants_decompose_mesh.h (from rev 14301, seismo/3D/SPECFEM3D_SESAME/trunk/UTILS/external_mesh/pre_meshfem3D_asteroid_subdivide/constants_decompose_mesh.h)
===================================================================
--- seismo/3D/SPECFEM3D_SESAME/trunk/UTILS/external_mesh/decompose_mesh/constants_decompose_mesh.h (rev 0)
+++ seismo/3D/SPECFEM3D_SESAME/trunk/UTILS/external_mesh/decompose_mesh/constants_decompose_mesh.h 2009-03-12 17:42:51 UTC (rev 14302)
@@ -0,0 +1,11 @@
+! 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
Deleted: seismo/3D/SPECFEM3D_SESAME/trunk/UTILS/external_mesh/decompose_mesh/constants_pre_meshfem3D.h
===================================================================
--- seismo/3D/SPECFEM3D_SESAME/trunk/UTILS/external_mesh/pre_meshfem3D_asteroid_subdivide/constants_pre_meshfem3D.h 2009-03-12 16:31:24 UTC (rev 14300)
+++ seismo/3D/SPECFEM3D_SESAME/trunk/UTILS/external_mesh/decompose_mesh/constants_pre_meshfem3D.h 2009-03-12 17:42:51 UTC (rev 14302)
@@ -1,11 +0,0 @@
-! 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
Copied: seismo/3D/SPECFEM3D_SESAME/trunk/UTILS/external_mesh/decompose_mesh/decompose_mesh.f90 (from rev 14301, seismo/3D/SPECFEM3D_SESAME/trunk/UTILS/external_mesh/pre_meshfem3D_asteroid_subdivide/decompose_mesh.f90)
===================================================================
--- seismo/3D/SPECFEM3D_SESAME/trunk/UTILS/external_mesh/decompose_mesh/decompose_mesh.f90 (rev 0)
+++ seismo/3D/SPECFEM3D_SESAME/trunk/UTILS/external_mesh/decompose_mesh/decompose_mesh.f90 2009-03-12 17:42:51 UTC (rev 14302)
@@ -0,0 +1,252 @@
+
+ program decompose_mesh
+
+ use part_decompose_mesh
+ implicit none
+
+ include './constants_decompose_mesh.h'
+ include "./scotchf.h"
+
+ integer, parameter :: nparts = 8
+
+ 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 :: vwgt
+ integer, dimension(:), allocatable :: adjwgt
+ integer, dimension(5) :: metis_options
+
+ 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 :: nb_materials
+ double precision, dimension(:), allocatable :: cs
+ integer, dimension(:), allocatable :: num_material
+
+ integer(long) :: nsize ! Max number of elements that contain the same node.
+ integer :: edgecut
+ integer :: nb_edges
+
+ integer :: ispec, inode
+ integer :: wgtflag
+ integer :: num_start
+ 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
+ character(len=256) :: prname
+
+ logical, dimension(:), allocatable :: mask_nodes_elmnts
+ integer, dimension(:), allocatable :: used_nodes_elmnts
+
+!!!! NL NL for SCOTCH partitioner
+ 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
+!!!! NL NL
+
+
+ ngnod = esize
+
+! read the elements, material and nodes files
+ open(unit=98, file='../model_asteroid_subdivide/mesh', status='old', form='formatted')
+ read(98,*) nspec
+ allocate(elmnts(esize,nspec))
+ do ispec = 1, nspec
+ read(98,*) elmnts(1,ispec), elmnts(2,ispec), elmnts(3,ispec), elmnts(4,ispec), &
+ elmnts(5,ispec), elmnts(6,ispec), elmnts(7,ispec), elmnts(8,ispec)
+ end do
+ close(98)
+
+ open(unit=98, file='../model_asteroid_subdivide/mat', status='old', form='formatted')
+ allocate(mat(nspec))
+ do ispec = 1, nspec
+ read(98,*) mat(ispec)
+ end do
+ close(98)
+
+ open(unit=98, file='../model_asteroid_subdivide/nodes_coords', status='old', form='formatted')
+ read(98,*) nnodes
+ allocate(nodes_coords(3,nnodes))
+ do inode = 1, nnodes
+ read(98,*) nodes_coords(1,inode), nodes_coords(2,inode), nodes_coords(3,inode)
+ end do
+ close(98)
+
+ 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
+ nsize = maxval(used_nodes_elmnts(:))
+ sup_neighbour = ngnod * nsize - (ngnod + (ngnod/2 - 1)*nfaces)
+ print*, 'nsize = ',nsize, 'sup_neighbour = ', sup_neighbour
+
+ do inode = 1, nnodes
+ if (.not. mask_nodes_elmnts(inode)) then
+ stop 'ERROR : nodes not used.'
+ endif
+ enddo
+! if (maxval(used_nodes_elmnts(:))>nsize) then
+! stop 'ERROR : increase nsize or modify the mesh.'
+! endif
+
+ elmnts(:,:) = elmnts(:,:) - 1
+
+ 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*, 'max_neighbour = ',max_neighbour
+
+
+ ! elmnts(:,:) = elmnts(:,:) + 1
+ ! adjncy(:) = adjncy(:) + 1
+ ! xadj(:) = xadj(:) + 1
+! allocate(vwgt(0:nspec-1))
+ nb_edges = xadj(nspec+1)
+! allocate(adjwgt(0:nb_edges-1))
+! vwgt(:) = 1
+! adjwgt(:) = 1
+
+! metis_options(1) = 0
+! metis_options(2) = 3
+! metis_options(3) = 1
+! metis_options(4) = 1
+! metis_options(5) = 0
+
+
+! num_start = 0
+! wgtflag = 0
+
+ allocate(part(1:nspec))
+
+
+! Old metis partitioning
+! call METIS_PartGraphRecursive(nspec, xadj(1), adjncy(1), vwgt(0), adjwgt(0), wgtflag, num_start, nparts, &
+! metis_options, edgecut, part(1));
+
+! 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
+
+ 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
+
+
+! local number of each element for each partition
+ call Construct_glob2loc_elmnts(nspec, part, nparts, glob2loc_elmnts)
+
+! local number of each node for each partition
+ call Construct_glob2loc_nodes(nspec, nnodes,nsize, nnodes_elmnts, nodes_elmnts, part, nparts, &
+ glob2loc_nodes_nparts, glob2loc_nodes_parts, glob2loc_nodes)
+
+ nb_materials = 1
+ allocate(cs(nb_materials))
+ allocate(num_material(nspec))
+ cs(:) = 1000.d0
+ num_material(:) = 1
+
+ call Construct_interfaces(nspec, nparts, sup_neighbour, part, elmnts, xadj, adjncy, tab_interfaces, &
+ tab_size_interfaces, ninterfaces, nb_materials, cs, num_material)
+
+ allocate(my_interfaces(0:ninterfaces-1))
+ allocate(my_nb_interfaces(0:ninterfaces-1))
+
+
+
+ do ipart = 0, nparts-1
+
+ !write(prname, "('/Database',i5.5)") ipart
+ write(prname, "(i6.6,'_Database')") ipart
+ open(unit=15,file='./OUTPUT_FILES/proc'//prname,status='unknown', action='write', form='formatted')
+
+ call write_glob2loc_nodes_database(15, ipart, nnodes_loc, nodes_coords, glob2loc_nodes_nparts, glob2loc_nodes_parts, &
+ glob2loc_nodes, nnodes, 1)
+ call write_partition_database(15, ipart, nspec_loc, nspec, elmnts, glob2loc_elmnts, glob2loc_nodes_nparts, &
+ glob2loc_nodes_parts, glob2loc_nodes, part, num_material, ngnod, 1)
+
+ 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)
+ 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, num_material, ngnod, 2)
+
+ call Write_interfaces_database(15, tab_interfaces, tab_size_interfaces, nparts, ipart, ninterfaces, &
+ my_ninterface, my_interfaces, my_nb_interfaces, glob2loc_elmnts, glob2loc_nodes_nparts, glob2loc_nodes_parts, &
+ glob2loc_nodes, 1)
+ write(15,*) my_ninterface, maxval(my_nb_interfaces)
+ call Write_interfaces_database(15, tab_interfaces, tab_size_interfaces, nparts, ipart, ninterfaces, &
+ my_ninterface, my_interfaces, my_nb_interfaces, glob2loc_elmnts, glob2loc_nodes_nparts, glob2loc_nodes_parts, &
+ glob2loc_nodes, 2)
+
+ close(15)
+
+ enddo
+
+ end program decompose_mesh
+
Copied: seismo/3D/SPECFEM3D_SESAME/trunk/UTILS/external_mesh/decompose_mesh/part_decompose_mesh.f90 (from rev 14301, seismo/3D/SPECFEM3D_SESAME/trunk/UTILS/external_mesh/pre_meshfem3D_asteroid_subdivide/part_decompose_mesh.f90)
===================================================================
--- seismo/3D/SPECFEM3D_SESAME/trunk/UTILS/external_mesh/decompose_mesh/part_decompose_mesh.f90 (rev 0)
+++ seismo/3D/SPECFEM3D_SESAME/trunk/UTILS/external_mesh/decompose_mesh/part_decompose_mesh.f90 2009-03-12 17:42:51 UTC (rev 14302)
@@ -0,0 +1,622 @@
+module part_decompose_mesh
+
+ implicit none
+
+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)
+
+ include './constants_decompose_mesh.h'
+
+ 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
+
+ integer :: i, j, k, l, m, nb_edges
+ logical :: is_neighbour
+ integer :: num_node, n
+ integer :: elem_base, elem_target
+ integer :: connectivity
+
+
+ !allocate(xadj(0:nelmnts))
+ xadj(:) = 0
+ !allocate(adjncy(0:max_neighbour*nelmnts-1))
+ adjncy(:) = 0
+ !allocate(nnodes_elmnts(0:nnodes-1))
+ nnodes_elmnts(:) = 0
+ !allocate(nodes_elmnts(0:nsize*nnodes-1))
+ 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
+
+ !print *, 'nnodes_elmnts'
+
+ ! 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, nparts, glob2loc_elmnts)
+
+ include './constants_decompose_mesh.h'
+
+ integer(long), intent(in) :: nelmnts
+ integer, intent(in) :: nparts
+ integer, dimension(0:nelmnts-1), intent(in) :: part
+ integer, dimension(:), pointer :: glob2loc_elmnts
+
+ integer :: num_glob, num_part
+ integer, dimension(0:nparts-1) :: num_loc
+
+
+ allocate(glob2loc_elmnts(0:nelmnts-1))
+
+ do num_part = 0, nparts-1
+ num_loc(num_part) = 0
+
+ end do
+
+ do num_glob = 0, nelmnts-1
+ num_part = part(num_glob)
+ 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, nparts, &
+ glob2loc_nodes_nparts, glob2loc_nodes_parts, glob2loc_nodes)
+
+ include './constants_decompose_mesh.h'
+
+ integer(long), intent(in) :: nelmnts, nsize
+ integer, intent(in) :: nnodes, nparts
+ 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
+ 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.
+ ! No interface between acoustic and elastic elements.
+ !--------------------------------------------------
+ subroutine Construct_interfaces(nelmnts, nparts, sup_neighbour, part, elmnts, xadj, adjncy, &
+ tab_interfaces, tab_size_interfaces, ninterfaces, nb_materials, cs_material, num_material)
+
+ include './constants_decompose_mesh.h'
+
+ integer, intent(in) :: 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
+ double precision, dimension(1:nb_materials), intent(in) :: cs_material
+ integer, intent(in) :: nb_materials
+
+
+ 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
+
+ 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
+
+ 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 ( cs_material(num_material(el+1)) < TINYVAL) then
+ is_acoustic_el = .true.
+ else
+ is_acoustic_el = .false.
+ end if
+ do el_adj = xadj(el), xadj(el+1)-1
+ if ( cs_material(num_material(adjncy(el_adj)+1)) < TINYVAL) then
+ is_acoustic_el_adj = .true.
+ 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
+ num_edge = num_edge + 1
+
+ end if
+ end do
+ end if
+ end do
+ 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
+
+ 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 ( cs_material(num_material(el+1)) < TINYVAL) then
+ is_acoustic_el = .true.
+ else
+ is_acoustic_el = .false.
+ end if
+ do el_adj = xadj(el), xadj(el+1)-1
+ if ( cs_material(num_material(adjncy(el_adj)+1)) < TINYVAL) then
+ is_acoustic_el_adj = .true.
+ 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
+
+
+
+ !--------------------------------------------------
+ ! 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)
+
+ integer, intent(in) :: IIN_database
+ integer, intent(in) :: nnodes, iproc, num_phase
+ integer, intent(inout) :: npgeo
+
+ 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
+ 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
+ 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 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.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(:) :: 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
+ nspec = 0
+
+ do i = 0, nelmnts-1
+ if ( part(i) == iproc ) then
+ nspec = nspec + 1
+
+ end if
+ end do
+
+ else
+ 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
+ write(IIN_database,*) glob2loc_elmnts(i)+1, num_modele(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, nparts, iproc, ninterfaces, &
+ my_ninterface, my_interfaces, my_nb_interfaces, glob2loc_elmnts, glob2loc_nodes_nparts, glob2loc_nodes_parts, &
+ glob2loc_nodes, num_phase)
+
+ integer, intent(in) :: IIN_database
+ integer, intent(in) :: iproc
+ integer, intent(in) :: nparts
+ integer, intent(in) :: ninterfaces
+ 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
+
+ num_interface = 0
+
+ if ( num_phase == 1 ) then
+
+ my_interfaces(:) = 0
+ my_nb_interfaces(:) = 0
+
+ do i = 0, nparts-1
+ do j = i+1, nparts-1
+ if ( (tab_size_interfaces(num_interface) < tab_size_interfaces(num_interface+1)) .and. &
+ (i == iproc .or. j == iproc) ) then
+ my_interfaces(num_interface) = 1
+ 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
+
+ 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
+
+ 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)
+ 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)
+ 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)
+ 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
+
+ end if
+
+ num_interface = num_interface + 1
+ end do
+ end do
+
+ end if
+
+ end subroutine write_interfaces_database
+
+end module part_decompose_mesh
+
Deleted: seismo/3D/SPECFEM3D_SESAME/trunk/UTILS/external_mesh/decompose_mesh/part_pre_meshfem3D.f90
===================================================================
--- seismo/3D/SPECFEM3D_SESAME/trunk/UTILS/external_mesh/pre_meshfem3D_asteroid_subdivide/part_pre_meshfem3D.f90 2009-03-12 16:31:24 UTC (rev 14300)
+++ seismo/3D/SPECFEM3D_SESAME/trunk/UTILS/external_mesh/decompose_mesh/part_pre_meshfem3D.f90 2009-03-12 17:42:51 UTC (rev 14302)
@@ -1,622 +0,0 @@
-module part_pre_meshfem3D
-
- implicit none
-
-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)
-
- include './constants_pre_meshfem3D.h'
-
- 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
-
- integer :: i, j, k, l, m, nb_edges
- logical :: is_neighbour
- integer :: num_node, n
- integer :: elem_base, elem_target
- integer :: connectivity
-
-
- !allocate(xadj(0:nelmnts))
- xadj(:) = 0
- !allocate(adjncy(0:max_neighbour*nelmnts-1))
- adjncy(:) = 0
- !allocate(nnodes_elmnts(0:nnodes-1))
- nnodes_elmnts(:) = 0
- !allocate(nodes_elmnts(0:nsize*nnodes-1))
- 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
-
- !print *, 'nnodes_elmnts'
-
- ! 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, nparts, glob2loc_elmnts)
-
- include './constants_pre_meshfem3D.h'
-
- integer(long), intent(in) :: nelmnts
- integer, intent(in) :: nparts
- integer, dimension(0:nelmnts-1), intent(in) :: part
- integer, dimension(:), pointer :: glob2loc_elmnts
-
- integer :: num_glob, num_part
- integer, dimension(0:nparts-1) :: num_loc
-
-
- allocate(glob2loc_elmnts(0:nelmnts-1))
-
- do num_part = 0, nparts-1
- num_loc(num_part) = 0
-
- end do
-
- do num_glob = 0, nelmnts-1
- num_part = part(num_glob)
- 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, nparts, &
- glob2loc_nodes_nparts, glob2loc_nodes_parts, glob2loc_nodes)
-
- include './constants_pre_meshfem3D.h'
-
- integer(long), intent(in) :: nelmnts, nsize
- integer, intent(in) :: nnodes, nparts
- 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
- 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.
- ! No interface between acoustic and elastic elements.
- !--------------------------------------------------
- subroutine Construct_interfaces(nelmnts, nparts, sup_neighbour, part, elmnts, xadj, adjncy, &
- tab_interfaces, tab_size_interfaces, ninterfaces, nb_materials, cs_material, num_material)
-
- include './constants_pre_meshfem3D.h'
-
- integer, intent(in) :: 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
- double precision, dimension(1:nb_materials), intent(in) :: cs_material
- integer, intent(in) :: nb_materials
-
-
- 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
-
- 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
-
- 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 ( cs_material(num_material(el+1)) < TINYVAL) then
- is_acoustic_el = .true.
- else
- is_acoustic_el = .false.
- end if
- do el_adj = xadj(el), xadj(el+1)-1
- if ( cs_material(num_material(adjncy(el_adj)+1)) < TINYVAL) then
- is_acoustic_el_adj = .true.
- 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
- num_edge = num_edge + 1
-
- end if
- end do
- end if
- end do
- 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
-
- 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 ( cs_material(num_material(el+1)) < TINYVAL) then
- is_acoustic_el = .true.
- else
- is_acoustic_el = .false.
- end if
- do el_adj = xadj(el), xadj(el+1)-1
- if ( cs_material(num_material(adjncy(el_adj)+1)) < TINYVAL) then
- is_acoustic_el_adj = .true.
- 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
-
-
-
- !--------------------------------------------------
- ! 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)
-
- integer, intent(in) :: IIN_database
- integer, intent(in) :: nnodes, iproc, num_phase
- integer, intent(inout) :: npgeo
-
- 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
- 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
- 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 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_pre_meshfem3D.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(:) :: 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
- nspec = 0
-
- do i = 0, nelmnts-1
- if ( part(i) == iproc ) then
- nspec = nspec + 1
-
- end if
- end do
-
- else
- 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
- write(IIN_database,*) glob2loc_elmnts(i)+1, num_modele(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, nparts, iproc, ninterfaces, &
- my_ninterface, my_interfaces, my_nb_interfaces, glob2loc_elmnts, glob2loc_nodes_nparts, glob2loc_nodes_parts, &
- glob2loc_nodes, num_phase)
-
- integer, intent(in) :: IIN_database
- integer, intent(in) :: iproc
- integer, intent(in) :: nparts
- integer, intent(in) :: ninterfaces
- 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
-
- num_interface = 0
-
- if ( num_phase == 1 ) then
-
- my_interfaces(:) = 0
- my_nb_interfaces(:) = 0
-
- do i = 0, nparts-1
- do j = i+1, nparts-1
- if ( (tab_size_interfaces(num_interface) < tab_size_interfaces(num_interface+1)) .and. &
- (i == iproc .or. j == iproc) ) then
- my_interfaces(num_interface) = 1
- 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
-
- 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
-
- 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)
- 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)
- 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)
- 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
-
- end if
-
- num_interface = num_interface + 1
- end do
- end do
-
- end if
-
- end subroutine write_interfaces_database
-
-
-end module part_pre_meshfem3D
Deleted: seismo/3D/SPECFEM3D_SESAME/trunk/UTILS/external_mesh/decompose_mesh/pre_meshfem3D.f90
===================================================================
--- seismo/3D/SPECFEM3D_SESAME/trunk/UTILS/external_mesh/pre_meshfem3D_asteroid_subdivide/pre_meshfem3D.f90 2009-03-12 16:31:24 UTC (rev 14300)
+++ seismo/3D/SPECFEM3D_SESAME/trunk/UTILS/external_mesh/decompose_mesh/pre_meshfem3D.f90 2009-03-12 17:42:51 UTC (rev 14302)
@@ -1,258 +0,0 @@
-program pre_meshfem3D
-
- use part_pre_meshfem3D
- implicit none
-
- include './constants_pre_meshfem3D.h'
- include "./scotchf.h"
-
-
- integer, parameter :: nparts=8
-
-
- 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 :: vwgt
- integer, dimension(:), allocatable :: adjwgt
- integer, dimension(5) :: metis_options
-
- 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 :: nb_materials
- double precision, dimension(:), allocatable :: cs
- integer, dimension(:), allocatable :: num_material
-
- integer(long) :: nsize ! Max number of elements that contain the same node.
- integer :: edgecut
- integer :: nb_edges
-
- integer :: ispec, inode
- integer :: wgtflag
- integer :: num_start
- 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
- character(len=256) :: prname
-
- logical, dimension(:), allocatable :: mask_nodes_elmnts
- integer, dimension(:), allocatable :: used_nodes_elmnts
-
-!!!! NL NL for SCOTCH partitioner
- 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
-!!!! NL NL
-
-
- ngnod = esize
-
-! read the elements, material and nodes files
- open(unit=98, file='../model_asteroid_subdivide/mesh', status='old', form='formatted')
- read(98,*) nspec
- allocate(elmnts(esize,nspec))
- do ispec = 1, nspec
- read(98,*) elmnts(1,ispec), elmnts(2,ispec), elmnts(3,ispec), elmnts(4,ispec), &
- elmnts(5,ispec), elmnts(6,ispec), elmnts(7,ispec), elmnts(8,ispec)
- end do
- close(98)
-
- open(unit=98, file='../model_asteroid_subdivide/mat', status='old', form='formatted')
- allocate(mat(nspec))
- do ispec = 1, nspec
- read(98,*) mat(ispec)
- end do
- close(98)
-
- open(unit=98, file='../model_asteroid_subdivide/nodes_coords', status='old', form='formatted')
- read(98,*) nnodes
- allocate(nodes_coords(3,nnodes))
- do inode = 1, nnodes
- read(98,*) nodes_coords(1,inode), nodes_coords(2,inode), nodes_coords(3,inode)
- end do
- close(98)
-
- 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
- nsize = maxval(used_nodes_elmnts(:))
- sup_neighbour = ngnod * nsize - (ngnod + (ngnod/2 - 1)*nfaces)
- print*, 'nsize = ',nsize, 'sup_neighbour = ', sup_neighbour
-
- do inode = 1, nnodes
- if (.not. mask_nodes_elmnts(inode)) then
- stop 'ERROR : nodes not used.'
- endif
- enddo
-! if (maxval(used_nodes_elmnts(:))>nsize) then
-! stop 'ERROR : increase nsize or modify the mesh.'
-! endif
-
- elmnts(:,:) = elmnts(:,:) - 1
-
- 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*, 'max_neighbour = ',max_neighbour
-
-
- ! elmnts(:,:) = elmnts(:,:) + 1
- ! adjncy(:) = adjncy(:) + 1
- ! xadj(:) = xadj(:) + 1
-! allocate(vwgt(0:nspec-1))
- nb_edges = xadj(nspec+1)
-! allocate(adjwgt(0:nb_edges-1))
-! vwgt(:) = 1
-! adjwgt(:) = 1
-
-! metis_options(1) = 0
-! metis_options(2) = 3
-! metis_options(3) = 1
-! metis_options(4) = 1
-! metis_options(5) = 0
-
-
-! num_start = 0
-! wgtflag = 0
-
- allocate(part(1:nspec))
-
-
-! Old metis partitioning
-! call METIS_PartGraphRecursive(nspec, xadj(1), adjncy(1), vwgt(0), adjwgt(0), wgtflag, num_start, nparts, &
-! metis_options, edgecut, part(1));
-
-! 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
-
- 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
-
-
-! local number of each element for each partition
- call Construct_glob2loc_elmnts(nspec, part, nparts, glob2loc_elmnts)
-
-! local number of each node for each partition
- call Construct_glob2loc_nodes(nspec, nnodes,nsize, nnodes_elmnts, nodes_elmnts, part, nparts, &
- glob2loc_nodes_nparts, glob2loc_nodes_parts, glob2loc_nodes)
-
- nb_materials = 1
- allocate(cs(nb_materials))
- allocate(num_material(nspec))
- cs(:) = 1000.d0
- num_material(:) = 1
-
- call Construct_interfaces(nspec, nparts, sup_neighbour, part, elmnts, xadj, adjncy, tab_interfaces, &
- tab_size_interfaces, ninterfaces, nb_materials, cs, num_material)
-
- allocate(my_interfaces(0:ninterfaces-1))
- allocate(my_nb_interfaces(0:ninterfaces-1))
-
-
-
- do ipart = 0, nparts-1
-
- !write(prname, "('/Database',i5.5)") ipart
- write(prname, "(i6.6,'_Database')") ipart
- open(unit=15,file='./OUTPUT_FILES/proc'//prname,status='unknown', action='write', form='formatted')
-
- call write_glob2loc_nodes_database(15, ipart, nnodes_loc, nodes_coords, glob2loc_nodes_nparts, glob2loc_nodes_parts, &
- glob2loc_nodes, nnodes, 1)
- call write_partition_database(15, ipart, nspec_loc, nspec, elmnts, glob2loc_elmnts, glob2loc_nodes_nparts, &
- glob2loc_nodes_parts, glob2loc_nodes, part, num_material, ngnod, 1)
-
- 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)
- 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, num_material, ngnod, 2)
-
- call Write_interfaces_database(15, tab_interfaces, tab_size_interfaces, nparts, ipart, ninterfaces, &
- my_ninterface, my_interfaces, my_nb_interfaces, glob2loc_elmnts, glob2loc_nodes_nparts, glob2loc_nodes_parts, &
- glob2loc_nodes, 1)
- write(15,*) my_ninterface, maxval(my_nb_interfaces)
- call Write_interfaces_database(15, tab_interfaces, tab_size_interfaces, nparts, ipart, ninterfaces, &
- my_ninterface, my_interfaces, my_nb_interfaces, glob2loc_elmnts, glob2loc_nodes_nparts, glob2loc_nodes_parts, &
- glob2loc_nodes, 2)
-
-
- close(15)
-
- end do
-
-
-end program pre_meshfem3D
-
-
-
-
Copied: seismo/3D/SPECFEM3D_SESAME/trunk/UTILS/external_mesh/files_needed_asteroid (from rev 14300, seismo/3D/SPECFEM3D_SESAME/trunk/UTILS/external_mesh/files_needed_for_asteroid_simulation)
Property changes on: seismo/3D/SPECFEM3D_SESAME/trunk/UTILS/external_mesh/files_needed_asteroid
___________________________________________________________________
Name: svn:mergeinfo
+
Deleted: seismo/3D/SPECFEM3D_SESAME/trunk/UTILS/external_mesh/readme
===================================================================
--- seismo/3D/SPECFEM3D_SESAME/trunk/UTILS/external_mesh/readme 2009-03-12 17:37:00 UTC (rev 14301)
+++ seismo/3D/SPECFEM3D_SESAME/trunk/UTILS/external_mesh/readme 2009-03-12 17:42:51 UTC (rev 14302)
@@ -1,53 +0,0 @@
-# A run using an external mesh is divided in five steps.
-#
-#
-# 1- we provide three files describing the model. There is an example of those in model_asteroid_subdivide.
-# - nodes_coords : this file contains the number of nodes on the first line, followed by the (x,y,z) coordinates for each node in meters.
-# - mesh : this file contains the number of elements on the first line, followed by the list of nodes (8 per element) that compose those elements. The ordering is the same as in specfem3D.
-# - mat : this file associates an element with its material (an integer).
-#
-# Those files were based on HOMO_3D_lisse_300, which is the abaqus export of a CUBIT mesh courtesy of Celine Blitz (celine.blitz<AT>univ-pau<DOT>fr) representing the EROS asteroid, and extracted using tail, head, and awk commands.
-#
-# A code named sub.f90 can be used to subdivide the mesh. It should ultimately be included in meshfem3D (see TODO list), as if the decimation is done prior to partitioning, METIS will have a hard time partioning it (NSUB in file constants.h is the number by which we decimate in all deirections).
-# Be careful with constants provided in constants_pre_meshfem3D.h, nsize and max_neighbour. There is currently no check as to whether the given mesh will be within those bounds. We should adress this issue shortly.
-#
-#
-# 2- we partition this mesh using the METIS software. pre_meshfem_asteroid_subdivide contains a stand alone program that makes the call to METIS, and extracts the faces/edges/nodes which NGLL points will fill the communication buffers. This should ultimately be done inside meshfem3D (see TODO list).
-# Compile with :
-# $ ./compil_all_sh
-# Run with :
-# $ ./a.out
-#
-# It writes files in OUTPUT_FILES named after the rank of the MPI process that will process this part. Those files must be copied into the OUTPUT_FILES where xmeshfem3D will be run.
-# In pre_meshfem3D.f90, the number of partitions should be indicated (nparts), and the path where resides the mesh, nodes_coords and mat files.
-#
-#
-# 3- we compile and run xmeshfem3D. The (rho,vp,vs) for the materials are hard written in create_region_mesh.f90, and should be modified accordingly. This would be better if it were read from a file (see TODO list). Prior to compilation, a close attention should be paid to the constants.h file. Several parameters are set inside this file, while they should better be read from a Par_file (see TODO list). To get to those parameters, look for string "nlegoff". There is a description for each of those parameters. A constants.h file is provided with parameters set for the asteroid simulation.
-#
-# Compile with :
-# $ ./configure FC= MPIFC=
-# $ make
-# Run with :
-# $ mpirun -n 96 ./xmeshfem3D
-#
-#
-# 4- we compile and run xspecfem3D. Again, A DATA/CMTSOLUTION file, a DATA/STATION and a DATA/PAr_file are provided for the asteroid run, while the format of those files should be changed.
-# - CMTSOLUTION : only four parameters are used : half duration is the main frequency (we run Ricker sources at the moment), latitude is y, longitude is x, depth is z (all three in meters).
-# - STATIONS : for each station, there is the station name, the network, y, x, a dummy value and z (x,y,z in meters).
-# - Par_file : practically no parameters are read from that file, except for NTSTEP_BETWEEN_FRAMES, USE_HIGHRES_FOR_MOVIES, LOCAL_PATH, NTSTEP_BETWEEN_OUTPUT_INFO. Others must be put to false in case of a run using an external mesh. It should really changed, and a new Par_file should gather info from the od Par_file along with the parametrs defined in constants.h.
-#
-# Compile with :
-# $ make xspecfem3D
-# Run with :
-# $ mpirun -n 96 ./xspecfem3D
-#
-#
-# 5- We create an OpenDX movie/shakemap using ./xcreate_movie_AVS_DX. Some parameters must be specified inside the file (prior to compilation).
-# - POWER_SCALING should be fixed to whatever suits the user. 0.13 was a good value for the asteroid shakemap, but 0.20 looked better for the movie.
-# - MUTE_SOURCE is used because we realized that the asteroid shakemap looked really bad due to the higher than average values near the source.
-# - RADIUS_TO_MUTE is the distance from the source below which the shakemap is nullified.
-# - X/Y/Z_SOURCE_EXT_MESH are the coordinates of the source.
-# - NSPEC_SURFACE_EXT_MESH is the number faces of elements that compose the envelope. Is printed as an output of specfem3D.
-#
-# Again, those parameters should be read from a PAr_file.
-#
Copied: seismo/3D/SPECFEM3D_SESAME/trunk/UTILS/external_mesh/readme_external_mesh (from rev 14300, seismo/3D/SPECFEM3D_SESAME/trunk/UTILS/external_mesh/readme)
===================================================================
--- seismo/3D/SPECFEM3D_SESAME/trunk/UTILS/external_mesh/readme_external_mesh (rev 0)
+++ seismo/3D/SPECFEM3D_SESAME/trunk/UTILS/external_mesh/readme_external_mesh 2009-03-12 17:42:51 UTC (rev 14302)
@@ -0,0 +1,53 @@
+# A run using an external mesh is divided in five steps.
+#
+#
+# 1- we provide three files describing the model. There is an example of those in model_asteroid_subdivide.
+# - nodes_coords : this file contains the number of nodes on the first line, followed by the (x,y,z) coordinates for each node in meters.
+# - mesh : this file contains the number of elements on the first line, followed by the list of nodes (8 per element) that compose those elements. The ordering is the same as in specfem3D.
+# - mat : this file associates an element with its material (an integer).
+#
+# Those files were based on HOMO_3D_lisse_300, which is the abaqus export of a CUBIT mesh courtesy of Celine Blitz (celine.blitz<AT>univ-pau<DOT>fr) representing the EROS asteroid, and extracted using tail, head, and awk commands.
+#
+# A code named sub.f90 can be used to subdivide the mesh. It should ultimately be included in meshfem3D (see TODO list), as if the decimation is done prior to partitioning, METIS will have a hard time partioning it (NSUB in file constants.h is the number by which we decimate in all deirections).
+# Be careful with constants provided in constants_pre_meshfem3D.h, nsize and max_neighbour. There is currently no check as to whether the given mesh will be within those bounds. We should adress this issue shortly.
+#
+#
+# 2- we partition this mesh using the METIS software. pre_meshfem_asteroid_subdivide contains a stand alone program that makes the call to METIS, and extracts the faces/edges/nodes which NGLL points will fill the communication buffers. This should ultimately be done inside meshfem3D (see TODO list).
+# Compile with :
+# $ ./compil_all_sh
+# Run with :
+# $ ./a.out
+#
+# It writes files in OUTPUT_FILES named after the rank of the MPI process that will process this part. Those files must be copied into the OUTPUT_FILES where xmeshfem3D will be run.
+# In pre_meshfem3D.f90, the number of partitions should be indicated (nparts), and the path where resides the mesh, nodes_coords and mat files.
+#
+#
+# 3- we compile and run xmeshfem3D. The (rho,vp,vs) for the materials are hard written in create_region_mesh.f90, and should be modified accordingly. This would be better if it were read from a file (see TODO list). Prior to compilation, a close attention should be paid to the constants.h file. Several parameters are set inside this file, while they should better be read from a Par_file (see TODO list). To get to those parameters, look for string "nlegoff". There is a description for each of those parameters. A constants.h file is provided with parameters set for the asteroid simulation.
+#
+# Compile with :
+# $ ./configure FC= MPIFC=
+# $ make
+# Run with :
+# $ mpirun -n 96 ./xmeshfem3D
+#
+#
+# 4- we compile and run xspecfem3D. Again, A DATA/CMTSOLUTION file, a DATA/STATION and a DATA/PAr_file are provided for the asteroid run, while the format of those files should be changed.
+# - CMTSOLUTION : only four parameters are used : half duration is the main frequency (we run Ricker sources at the moment), latitude is y, longitude is x, depth is z (all three in meters).
+# - STATIONS : for each station, there is the station name, the network, y, x, a dummy value and z (x,y,z in meters).
+# - Par_file : practically no parameters are read from that file, except for NTSTEP_BETWEEN_FRAMES, USE_HIGHRES_FOR_MOVIES, LOCAL_PATH, NTSTEP_BETWEEN_OUTPUT_INFO. Others must be put to false in case of a run using an external mesh. It should really changed, and a new Par_file should gather info from the od Par_file along with the parametrs defined in constants.h.
+#
+# Compile with :
+# $ make xspecfem3D
+# Run with :
+# $ mpirun -n 96 ./xspecfem3D
+#
+#
+# 5- We create an OpenDX movie/shakemap using ./xcreate_movie_AVS_DX. Some parameters must be specified inside the file (prior to compilation).
+# - POWER_SCALING should be fixed to whatever suits the user. 0.13 was a good value for the asteroid shakemap, but 0.20 looked better for the movie.
+# - MUTE_SOURCE is used because we realized that the asteroid shakemap looked really bad due to the higher than average values near the source.
+# - RADIUS_TO_MUTE is the distance from the source below which the shakemap is nullified.
+# - X/Y/Z_SOURCE_EXT_MESH are the coordinates of the source.
+# - NSPEC_SURFACE_EXT_MESH is the number faces of elements that compose the envelope. Is printed as an output of specfem3D.
+#
+# Again, those parameters should be read from a PAr_file.
+#
Property changes on: seismo/3D/SPECFEM3D_SESAME/trunk/UTILS/external_mesh/readme_external_mesh
___________________________________________________________________
Name: svn:mergeinfo
+
Copied: seismo/3D/SPECFEM3D_SESAME/trunk/analyze_CUBIT_Abaqus_mesh (from rev 14300, seismo/3D/SPECFEM3D_SESAME/trunk/UTILS/external_mesh/analyze_CUBIT_Abaqus_mesh)
Property changes on: seismo/3D/SPECFEM3D_SESAME/trunk/analyze_CUBIT_Abaqus_mesh
___________________________________________________________________
Name: svn:mergeinfo
+
More information about the CIG-COMMITS
mailing list