[cig-commits] r18359 - in seismo/3D/FAULT_SOURCE/branches/new_fault_db: decompose_mesh_SCOTCH src
ampuero at geodynamics.org
ampuero at geodynamics.org
Thu May 12 22:50:14 PDT 2011
Author: ampuero
Date: 2011-05-12 22:50:14 -0700 (Thu, 12 May 2011)
New Revision: 18359
Modified:
seismo/3D/FAULT_SOURCE/branches/new_fault_db/decompose_mesh_SCOTCH/decompose_mesh_SCOTCH.f90
seismo/3D/FAULT_SOURCE/branches/new_fault_db/decompose_mesh_SCOTCH/fault_scotch.f90
seismo/3D/FAULT_SOURCE/branches/new_fault_db/src/create_regions_mesh.f90
seismo/3D/FAULT_SOURCE/branches/new_fault_db/src/fault_object.f90
Log:
closing fault in fault_object (instead of fault_scotch)
Modified: seismo/3D/FAULT_SOURCE/branches/new_fault_db/decompose_mesh_SCOTCH/decompose_mesh_SCOTCH.f90
===================================================================
--- seismo/3D/FAULT_SOURCE/branches/new_fault_db/decompose_mesh_SCOTCH/decompose_mesh_SCOTCH.f90 2011-05-13 02:53:59 UTC (rev 18358)
+++ seismo/3D/FAULT_SOURCE/branches/new_fault_db/decompose_mesh_SCOTCH/decompose_mesh_SCOTCH.f90 2011-05-13 05:50:14 UTC (rev 18359)
@@ -389,7 +389,8 @@
if( nspec2D_moho > 0 ) print*, ' nspec2D_moho = ', nspec2D_moho
call read_fault_files(localpath_name)
- call close_faults(nodes_coords,elmnts,nspec,nnodes,esize)
+!jpa: keep the fault open, we'll close it in fault_object.f90
+! call close_faults(nodes_coords,elmnts,nspec,nnodes,esize)
end subroutine read_mesh_files
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-13 02:53:59 UTC (rev 18358)
+++ seismo/3D/FAULT_SOURCE/branches/new_fault_db/decompose_mesh_SCOTCH/fault_scotch.f90 2011-05-13 05:50:14 UTC (rev 18359)
@@ -157,25 +157,38 @@
do j = 1,size(ispec1)
do k1=1,esize
- iglob1 = elmnts(k1,ispec1(j))
- xyz_1 = nodes_coords(:,iglob1)
- xyz = xyz_2-xyz_1
- dist = xyz(1)*xyz(1) + xyz(2)*xyz(2) + xyz(3)*xyz(3)
- dist = sqrt(dist)
- !pgb : dist has to be > 0 in case we want to avoid re sticking nodes with
- !pgb : the same positions already , like nodes at fault edeges.
- !pgb : Or nodes that have been sticked already.
-
- if ((0.0d0 < dist) .and. (dist <= FAULT_GAP_TOLERANCE)) then
- xyz = (xyz_1 + xyz_2)*0.5d0
- nodes_coords(:,iglob2) = xyz
- nodes_coords(:,iglob1) = xyz
- found_it = .true.
- exit
- endif
+ iglob1 = elmnts(k1,ispec1(j))
+ xyz_1 = nodes_coords(:,iglob1)
+ xyz = xyz_2-xyz_1
+ dist = xyz(1)*xyz(1) + xyz(2)*xyz(2) + xyz(3)*xyz(3)
+ dist = sqrt(dist)
+
+ !pgb : dist has to be > 0 in case we want to avoid re sticking nodes with
+ !pgb : the same positions already , like nodes at fault edeges.
+ !pgb : Or nodes that have been sticked already.
+
+ !jpa: Closing nodes that are already closed is not a problem
+ !jpa: I process them again to leave the loop as early as possible
+ !jpa: and to test if a facing node was found (See below).
+
+ !if ((0.0d0 < dist) .and. (dist <= FAULT_GAP_TOLERANCE)) then
+ if (dist <= FAULT_GAP_TOLERANCE) then
+ xyz = (xyz_1 + xyz_2)*0.5d0
+ nodes_coords(:,iglob2) = xyz
+ nodes_coords(:,iglob1) = xyz
+ found_it = .true.
+ exit
+ endif
+ enddo
- enddo
- if (found_it) exit
+ if (found_it) then
+ exit
+ else
+ ! If the fault is very complicated (non-planar) the meshes of the two fault sides
+ ! can be very unstructured and might not match (unless you enforce it in CUBIT).
+ ! That is unlikely because the fault gap is tiny, but we never know.
+ stop 'Inconsistent fault mesh: corresponding node in the other fault face was not found'
+ endif
enddo
enddo
Modified: seismo/3D/FAULT_SOURCE/branches/new_fault_db/src/create_regions_mesh.f90
===================================================================
--- seismo/3D/FAULT_SOURCE/branches/new_fault_db/src/create_regions_mesh.f90 2011-05-13 02:53:59 UTC (rev 18358)
+++ seismo/3D/FAULT_SOURCE/branches/new_fault_db/src/create_regions_mesh.f90 2011-05-13 05:50:14 UTC (rev 18359)
@@ -177,7 +177,8 @@
write(IMAIN,*) ' ...setting up jacobian '
endif
-
+!jpa: jacobians should be computed after closgin the fault (after fault_setup)
+!jpa: can we tolerate this small error?
call crm_ext_setup_jacobian(myrank, &
xstore,ystore,zstore,nspec, &
nodes_coords_ext_mesh,nnodes_ext_mesh,&
Modified: seismo/3D/FAULT_SOURCE/branches/new_fault_db/src/fault_object.f90
===================================================================
--- seismo/3D/FAULT_SOURCE/branches/new_fault_db/src/fault_object.f90 2011-05-13 02:53:59 UTC (rev 18358)
+++ seismo/3D/FAULT_SOURCE/branches/new_fault_db/src/fault_object.f90 2011-05-13 05:50:14 UTC (rev 18359)
@@ -147,25 +147,19 @@
subroutine fault_setup(ibool,nnodes_ext_mesh,nodes_coords_ext_mesh, &
xstore,ystore,zstore,nspec,nglob,myrank)
- integer, intent(in) :: nspec ! number of spectral elements in each block
- double precision, dimension(NGLLX,NGLLY,NGLLZ,nspec), intent(in) :: xstore,ystore,zstore
+ integer, intent(in) :: nspec
+ double precision, dimension(NGLLX,NGLLY,NGLLZ,nspec), intent(inout) :: xstore,ystore,zstore
integer, dimension(NGLLX,NGLLY,NGLLZ,nspec), intent(in) :: ibool
integer, intent(in) :: myrank
integer, intent(in) :: nglob
integer, intent(in) :: nnodes_ext_mesh
double precision, dimension(NDIM,nnodes_ext_mesh),intent(in) :: nodes_coords_ext_mesh
-
integer :: iflt
logical :: fault_exists
if (.not. allocated(fault_db)) return
-! 1. Generate node indexing (ibool) from original coordinates. Fault nodes are split in CUBIT.
-! 2. Slightly shift the coordinates of nodes on side 1 of the fault
-
-! to do: what happens to nodes at the edges of a fault ??? : still thinking.
-
fault_exists = .false.
do iflt=1,size(fault_db)
if (fault_db(iflt)%nspec>0) then
@@ -192,6 +186,8 @@
call setup_Kelvin_Voigt_eta(fault_db(iflt),nspec)
+ ! close the fault (modifies *store_dummy) to compute jacobians and normals
+ call close_fault(fault_db(iflt)) !,xstore,ystore,zstore,nspec)
call setup_normal_jacobian(fault_db(iflt),ibool,nspec,nglob,myrank)
enddo
@@ -199,6 +195,59 @@
end subroutine fault_setup
!=============================================================================================================
+! We only close *store_dummy
+! Fortunately only *store_dummy is needed to compute jacobians and normals
+subroutine close_fault(fdb) !,xstore,ystore,zstore,nspec)
+
+ use create_regions_mesh_ext_par, only: xstore_dummy,ystore_dummy,zstore_dummy
+
+ type(fault_db_type), intent(inout) :: fdb
+! integer, intent(in) :: nspec
+! double precision, dimension(NGLLX,NGLLY,NGLLZ,nspec), intent(inout) :: xstore,ystore,zstore
+
+ integer :: i,K1,K2
+
+ do i=1,fdb%nglob
+ K1 = fdb%ibulk1(i)
+ K2 = fdb%ibulk2(i)
+ xstore_dummy(K1) = 0.5d0*( xstore_dummy(K1) + xstore_dummy(K2) )
+ xstore_dummy(K2) = xstore_dummy(K1)
+ ystore_dummy(K1) = 0.5d0*( ystore_dummy(K1) + ystore_dummy(K2) )
+ ystore_dummy(K2) = ystore_dummy(K1)
+ zstore_dummy(K1) = 0.5d0*( zstore_dummy(K1) + zstore_dummy(K2) )
+ zstore_dummy(K2) = zstore_dummy(K1)
+ enddo
+
+! Closing *store is not easy ...
+! the following would work only if ispec1 and ispec2 were always facing each other
+! and if igll corresponded to the same node in both element faces
+!
+! do e=1,fdb%nspec
+! e1 = fdb%ispec1(e)
+! e2 = fdb%ispec1(e)
+! igll = 0
+! do i=1,NGLLX
+! do j=1,NGLLX
+! igll = igll + 1
+! i1=fdb%ijk1(1,igll,e)
+! j1=fdb%ijk1(2,igll,e)
+! k1=fdb%ijk1(3,igll,e)
+! i2=fdb%ijk2(1,igll,e)
+! j2=fdb%ijk2(2,igll,e)
+! k2=fdb%ijk2(3,igll,e)
+! xstore(i1,j1,k1,e1) = 0.5d0*( xstore(i1,j1,k1,e1) + xstore(i2,j2,k2,e2) )
+! xstore(i2,j2,k2,e2) = xstore(i1,j1,k1,e1)
+! ystore(i1,j1,k1,e1) = 0.5d0*( ystore(i1,j1,k1,e1) + ystore(i2,j2,k2,e2) )
+! ystore(i2,j2,k2,e2) = ystore(i1,j1,k1,e1)
+! zstore(i1,j1,k1,e1) = 0.5d0*( zstore(i1,j1,k1,e1) + zstore(i2,j2,k2,e2) )
+! zstore(i2,j2,k2,e2) = zstore(i1,j1,k1,e1)
+! enddo
+! enddo
+! enddo
+
+subroutine close_fault
+
+!=============================================================================================================
! looks for i,j,k indices of GLL points on boundary face
! determines element face by given CUBIT corners
! sets face id of reference element associated with this face
@@ -285,7 +334,7 @@
end subroutine
!===============================================================================================================
-! The lexicographic oredering of node coordinates
+! The lexicographic ordering of node coordinates
! guarantees that the fault nodes are
! consistently ordered on both sides of the fault,
! such that the K-th node of side 1 is facing the K-th node of side 2
More information about the CIG-COMMITS
mailing list