[cig-commits] r18498 - seismo/3D/FAULT_SOURCE/branches/new_fault_db/decompose_mesh_SCOTCH

[ampuero at geodynamics.org]

Author: ampuero
Date: 2011-05-30 23:04:26 -0700 (Mon, 30 May 2011)
New Revision: 18498

Modified:
   seismo/3D/FAULT_SOURCE/branches/new_fault_db/decompose_mesh_SCOTCH/fault_scotch.f90
Log:
redefined fault zone layer in scotch, now based on fault%inodes# to keep it thin

Modified: seismo/3D/FAULT_SOURCE/branches/new_fault_db/decompose_mesh_SCOTCH/fault_scotch.f90
===================================================================
--- seismo/3D/FAULT_SOURCE/branches/new_fault_db/decompose_mesh_SCOTCH/fault_scotch.f90	2011-05-31 05:27:02 UTC (rev 18497)
+++ seismo/3D/FAULT_SOURCE/branches/new_fault_db/decompose_mesh_SCOTCH/fault_scotch.f90	2011-05-31 06:04:26 UTC (rev 18498)
@@ -16,8 +16,7 @@
  
   integer, parameter :: long = SELECTED_INT_KIND(18)
 
-  double precision, parameter :: FAULT_GAP_TOLERANCE = 1.0d0 ! JPA: are you sure 1 meter is small enough?
-                                                             ! PGB: for a simple test is fine .For SCEC lower values.
+  double precision, parameter :: FAULT_GAP_TOLERANCE = 0.01d0
 
   public :: read_fault_files, fault_collect_elements, close_faults, write_fault_database, &
             save_nodes_coords, nodes_coords_open, faults
@@ -253,15 +252,15 @@
   integer, dimension(0:nelmnts-1), intent(inout)    :: part
 
 !LOCAL VARIABLES :
-  integer, dimension(0:nelmnts)                  :: xadj
-  integer, dimension(0:sup_neighbour*nelmnts-1)  :: adjncy
+!  integer, dimension(0:nelmnts)                  :: xadj
+!  integer, dimension(0:sup_neighbour*nelmnts-1)  :: adjncy
   integer, dimension(0:nnodes-1)                 :: nnodes_elmnts
   integer, dimension(0:nsize*nnodes-1)           :: nodes_elmnts
-  integer  :: max_neighbour       
+!  integer  :: max_neighbour       
 
 !SHILDING 
   integer  :: i,j, ipart,nproc_null,nproc_null_final
-  integer  :: el, el_1, el_2, k1, k2
+  integer  :: el, el_1, el_2, k1, k2, k,e,iflt,inode
   logical  :: is_repartitioned
   integer, dimension(:), allocatable :: elem_proc_null
 
@@ -296,9 +295,9 @@
     end do
 
   endif
-  call mesh2dual_ncommonnodes_fault(nelmnts, nnodes, nsize, sup_neighbour, &
-                                elmnts, xadj, adjncy, nnodes_elmnts, &
-                                nodes_elmnts, max_neighbour, 4, esize)
+ ! call mesh2dual_ncommonnodes_fault(nelmnts, nnodes, nsize, sup_neighbour, &
+ !                               elmnts, xadj, adjncy, nnodes_elmnts, &
+ !                               nodes_elmnts, max_neighbour, 4, esize)
   
     ! coupled elements
     !  ---------------
@@ -314,21 +313,52 @@
     !                1    2   
                  
     ! Allocating elements with double shield layer
-  print *, "Fault shield double-layer :"
-  do el = 0, nelmnts-1
-    if ( is_on_fault(el+1) ) then
-      part(el) = 0
-      do k1 = xadj(el), xadj(el+1) - 1
-        el_1 = adjncy(k1) 
-        part(el_1) = 0
-        do k2 = xadj(el_1), xadj(el_1+1) - 1
-          el_2 = adjncy(k2) 
-          part(el_2) = 0
-        enddo
-      enddo
-    endif
-  enddo
+!  print *, "Fault shield double-layer :"
+!  do el = 0, nelmnts-1
+!    if ( is_on_fault(el+1) ) then
+!      part(el) = 0
+!      do k1 = xadj(el), xadj(el+1) - 1
+!        el_1 = adjncy(k1) 
+!        part(el_1) = 0
+!        do k2 = xadj(el_1), xadj(el_1+1) - 1
+!          el_2 = adjncy(k2) 
+!          part(el_2) = 0
+!        enddo
+!      enddo
+!    endif
+!  enddo
 
+ ! list of elements per node
+ !   nnodes_elmnts(i) = number of elements containing node #i (i=0:nnodes-1)     
+ !   nodes_elmnts(nsize*i:nsize*i+nnodes_elmnts(i)-1) = index of elements containing node #i  (k=0:nsize*nnodes-1)
+  nnodes_elmnts(:) = 0
+  nodes_elmnts(:) = 0
+  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
+
+ ! fault zone layer = the set of elements that contain at least one fault node
+  print *, "Fault zone layer :"
+  do iflt=1,size(faults)
+    do e=1,faults(iflt)%nspec
+      do k=1,4
+
+        inode = faults(iflt)%inodes1(k,e)-1  ! node index, starting at 0
+        k1 = nsize*inode
+        k2 = k1 + nnodes_elmts(inode) -1
+        part( nodes_elmnts(k1:k2) ) = 0
+
+        inode = faults(iflt)%inodes2(k,e)-1 
+        k1 = nsize*inode
+        k2 = k1 + nnodes_elmts(inode) -1
+        part( nodes_elmnts(k1:k2) ) = 0
+
+      end do
+    end do
+  end do
+! jpa: if the previous section works then clean up obsolete stuff
+
   nproc_null_final = count( part == 0 )
   print *, nproc_null_final