[cig-commits] r18322 - in seismo/3D/FAULT_SOURCE/branches/new_fault_db: . CUBIT decompose_mesh_SCOTCH src
percygalvez at geodynamics.org
percygalvez at geodynamics.org
Thu May 5 12:52:23 PDT 2011
Author: percygalvez
Date: 2011-05-05 12:52:23 -0700 (Thu, 05 May 2011)
New Revision: 18322
Added:
seismo/3D/FAULT_SOURCE/branches/new_fault_db/CUBIT/
seismo/3D/FAULT_SOURCE/branches/new_fault_db/CUBIT/crack_fault.py
Removed:
seismo/3D/FAULT_SOURCE/branches/new_fault_db/src/fault_ibool.f90
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/decompose_mesh_SCOTCH/part_decompose_mesh_SCOTCH.f90
seismo/3D/FAULT_SOURCE/branches/new_fault_db/src/fault_object.f90
Log:
fault_object correct and SCOTCH
Added: seismo/3D/FAULT_SOURCE/branches/new_fault_db/CUBIT/crack_fault.py
===================================================================
--- seismo/3D/FAULT_SOURCE/branches/new_fault_db/CUBIT/crack_fault.py (rev 0)
+++ seismo/3D/FAULT_SOURCE/branches/new_fault_db/CUBIT/crack_fault.py 2011-05-05 19:52:23 UTC (rev 18322)
@@ -0,0 +1,118 @@
+#!/usr/bin/env python
+
+import cubit
+import boundary_definition
+import cubit2specfem3d
+import math
+import os
+import sys
+cubit.cmd('reset')
+
+km = 1000
+z_surf=0*km
+#### initializing coordinates x,y,z
+x=[] # fault
+y=[]
+z=[]
+
+xbulk=[] # bulk
+ybulk=[]
+zbulk=[]
+
+xbulk.append(-21*km) #x1
+xbulk.append(21*km) #x2
+xbulk.append(21*km) #x3
+xbulk.append(-21*km) #x4
+
+ybulk.append(-16.2*km) #y1
+ybulk.append(-16.2*km) #y2
+ybulk.append(16.2*km) #y3
+ybulk.append(16.2*km) #y4
+
+zbulk=[z_surf]*4
+
+x.append(-16*km) #x5
+x.append(0*km) #x6
+x.append(16*km) #x7
+x.append(0*km) #x8
+
+y.append(0.0) #y5
+y.append(0.5*km) #y6
+y.append(0.0) #y7
+y.append(-0.5*km) #y8
+
+
+z=[z_surf]*4
+
+#################### bulk ###########################################
+for i in range(len(xbulk)):
+ vert="create vertex x "+str(xbulk[i])+" y "+str(ybulk[i])+" z "+str(zbulk[i])
+ cubit.cmd(vert)
+
+################ Loading fault points profile#############################
+for i in range(len(x)):
+ vert="create vertex x "+str(x[i])+" y "+str(y[i])+" z "+str(z[i])
+ cubit.cmd(vert)
+
+
+################ creating fault domains #################################
+bulk1="create curve vertex 1 2"
+bulk2="create curve vertex 2 3"
+bulk3="create curve vertex 3 4"
+bulk4="create curve vertex 4 1"
+
+fault_up="create curve spline vertex 5 6 7"
+fault_down="create curve spline vertex 5 8 7"
+
+cubit.cmd(bulk1)
+cubit.cmd(bulk2)
+cubit.cmd(bulk3)
+cubit.cmd(bulk4)
+
+cubit.cmd(fault_up)
+cubit.cmd(fault_down)
+
+surface="create surface curve 1 2 3 4 5 6"
+cubit.cmd(surface)
+
+cubit.cmd("sweep surface 1 vector 0 0 -1 distance "+str(21*km))
+cubit.cmd("sweep curve 5 vector 0 0 -1 distance "+str(21*km))
+cubit.cmd("sweep curve 6 vector 0 0 -1 distance "+str(21*km))
+
+#####################################################
+cubit.cmd("imprint all")
+cubit.cmd("merge all")
+cubit.cmd("mesh surface 1")
+cubit.cmd("mesh volume 1")
+cubit.cmd("unmerge surface 2 3")
+
+########### Elements nodes ###############
+fault_A_elements_up="sideset 200 surface 2"
+fault_A_elements_down="sideset 201 surface 3"
+cubit.cmd(fault_A_elements_up)
+cubit.cmd(fault_A_elements_down)
+
+
+#nsets=cubit.get_sideset_id_list()
+### SIDESETS ARE ELEMENTSSSSSSSSSSSS #####
+# To get nodes on fault elments used : cubit.get_connectivity(entity_type,entity_id)
+# To get faces on fault elments used : cubit.get_sub_elements(entity_type,entity_id,dimension)
+
+
+for nset in nsets:
+ el=cubit.get_sub_elements("hex",nset,2)
+ print el
+
+# name=cubit.get_exodus_entity_name("sideset",nset)
+# id_sideset=cubit.get_exodus_id("hex",nset)
+# print name
+# print id_sideset
+
+
+
+
+
+
+
+
+
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-05 07:41:11 UTC (rev 18321)
+++ seismo/3D/FAULT_SOURCE/branches/new_fault_db/decompose_mesh_SCOTCH/decompose_mesh_SCOTCH.f90 2011-05-05 19:52:23 UTC (rev 18322)
@@ -401,7 +401,6 @@
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-05 07:41:11 UTC (rev 18321)
+++ seismo/3D/FAULT_SOURCE/branches/new_fault_db/decompose_mesh_SCOTCH/fault_scotch.f90 2011-05-05 19:52:23 UTC (rev 18322)
@@ -5,7 +5,7 @@
type fault_type
private
- integer :: nspec
+ integer :: nspec,eta
integer, dimension(:), pointer :: ispec1, ispec2, iface1, iface2
end type fault_type
@@ -26,7 +26,7 @@
Subroutine read_fault_files
- integer :: nbfaults, iflt, ier
+ integer :: nbfaults=0, iflt, ier
open(101,file='../DATA/FAULT/Par_file_faults.in',status='old',action='read',iostat=ier)
if (ier==0) then
@@ -34,15 +34,15 @@
else
print *, 'Par_file.in not found: assume no faults'
endif
- close(101)
if (nbfaults>0) then
allocate(faults(nbfaults))
do iflt = 1 , nbfaults
- call read_single_fault_file(faults(iflt),iflt)
+ call read_single_fault_file(faults(iflt),iflt)
enddo
endif
+ close(101)
end subroutine read_fault_files
@@ -55,12 +55,12 @@
integer,intent(in) :: ifault
character(len=5) :: NTchar
- integer :: i,j
+ integer :: i,j,ier
write(NTchar,'(I5)') ifault
NTchar = adjustl(NTchar)
- open(101,file='../DATA/FAULT/fault_elements_'//NTCHAR//'.dat', &
+ open(101,file='../DATA/FAULT/fault_elements_'//NTCHAR//'.dat',&
status='old',action='read',iostat=ier)
if( ier == 0 ) then
@@ -88,13 +88,13 @@
subroutine close_faults(nodes_coords,elmnts,nelmnts,nnodes,esize)
integer ,intent(in) :: nnodes, esize, nelmnts
- integer, dimension(3,nnodes), intent(in) :: nodes_coords
+ real(kind=CUSTOM_REAL), dimension(3,nnodes), intent(inout) :: nodes_coords
integer, dimension(esize,nelmnts), intent(in) :: elmnts
integer :: iflt
do iflt=1,size(faults)
- call close_fault_single(faults(iflt)%ispec1,faults(iflt)%ispec2, &
+ call close_fault_single(faults(iflt)%ispec1,faults(iflt)%ispec2,&
elmnts,nodes_coords,nnodes,esize,nelmnts)
enddo
@@ -106,9 +106,9 @@
integer ,intent(in) :: nnodes, esize, nelmnts
integer, dimension(esize,nelmnts), intent(in) :: elmnts
integer , dimension(:), intent(in) :: ispec1,ispec2
- real(kind=CUSTOM_REAL),dimension(3,nnodes),target, intent(inout) :: nodes_coords
+ real(kind=CUSTOM_REAL),dimension(3,nnodes), intent(inout) :: nodes_coords
- real(kind=CUSTOM_REAL), dimension(3),pointer :: xyz_1 =>null(),xyz_2=>null()
+ real(kind=CUSTOM_REAL), dimension(3) :: xyz_1, xyz_2
real(kind=CUSTOM_REAL), dimension(3) :: xyz
real(kind=CUSTOM_REAL) :: dist
@@ -120,13 +120,13 @@
iglob2 = elmnts(k2,ispec2(i))
found_it = .false.
- xyz_2 => nodes_coords(:,iglob2)
+ xyz_2 = nodes_coords(:,iglob2)
do j = 1,size(ispec1)
do k1=1,esize
iglob1 = elmnts(k1,ispec1(j))
- xyz_1 => nodes_coords(:,iglob1)
+ xyz_1 = nodes_coords(:,iglob1)
xyz = xyz_2-xyz_1
dist = xyz(1)*xyz(1) + xyz(2)*xyz(2) + xyz(3)*xyz(3)
Modified: seismo/3D/FAULT_SOURCE/branches/new_fault_db/decompose_mesh_SCOTCH/part_decompose_mesh_SCOTCH.f90
===================================================================
--- seismo/3D/FAULT_SOURCE/branches/new_fault_db/decompose_mesh_SCOTCH/part_decompose_mesh_SCOTCH.f90 2011-05-05 07:41:11 UTC (rev 18321)
+++ seismo/3D/FAULT_SOURCE/branches/new_fault_db/decompose_mesh_SCOTCH/part_decompose_mesh_SCOTCH.f90 2011-05-05 19:52:23 UTC (rev 18322)
@@ -524,7 +524,6 @@
integer, intent(in) :: IIN_database
integer, intent(in) :: nnodes, iproc, num_phase
integer, intent(inout) :: npgeo
- integer, intent(in),dimension(:) :: ispec1,ispec2
double precision, dimension(3,nnodes) :: nodes_coords
integer, dimension(:), pointer :: glob2loc_nodes_nparts
Deleted: seismo/3D/FAULT_SOURCE/branches/new_fault_db/src/fault_ibool.f90
===================================================================
--- seismo/3D/FAULT_SOURCE/branches/new_fault_db/src/fault_ibool.f90 2011-05-05 07:41:11 UTC (rev 18321)
+++ seismo/3D/FAULT_SOURCE/branches/new_fault_db/src/fault_ibool.f90 2011-05-05 19:52:23 UTC (rev 18322)
@@ -1,11 +0,0 @@
-subroutine fault_ibool()
-
-! number of fault nodes : nfspec
-! initial pointers : loc
-! xp,yp,zp= coordinates of fault elements.
-! ije=ke ,
-! xp=xcoor(ke) ,
-! yp=ycoor(ke) ,
-! zp=zcoor(ke) ,
-
-call get_global(nspec,xp(ije),yp(ije),zp(ije),fault_ibool,loc,ifseg,nfault_ibool,npointot,UTM_X_MIN,UTM_X_MAX)
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-05 07:41:11 UTC (rev 18321)
+++ seismo/3D/FAULT_SOURCE/branches/new_fault_db/src/fault_object.f90 2011-05-05 19:52:23 UTC (rev 18322)
@@ -37,7 +37,7 @@
type fault_db_type
private
- integer :: tag1,tag2,nspec=0,nglob=0
+ integer :: nspec=0,nglob=0
real(kind=CUSTOM_REAL) :: eta
integer, dimension(:), pointer:: ispec1, ispec2, ibulk1, ibulk2, iface1, iface2
real(kind=CUSTOM_REAL), dimension(:), pointer :: xcoordbulk1,ycoordbulk1,zcoordbulk1,xcoordbulk2,ycoordbulk2,zcoordbulk2
@@ -87,6 +87,7 @@
read(100,*) nb
allocate(fault_db(nb))
do i=1,nb
+ fault_db(i)%eta
enddo
else
write(6,*) 'File Par_file_faults.in does not exist '
@@ -137,14 +138,13 @@
fault_exists = .false.
+ call loading_coords_and_setup(fault_db,nglob,nspec,prname,myrank)
+
do iflt=1,size(fault_db)
- call loading_coords_and_setup(fault_db(iflt),nglob,nspec,prname,myrank)
if (fault_db(iflt)%nspec>0) fault_exists = .true.
enddo
if (fault_exists) then
- call crm_ext_setup_indexing_fault(ibool, &
- xtemp,ytemp,ztemp,nspec,nglob,npointot)
!-------------- Kelvin voigt damping -------------------------
KELVIN_VOIGT_DAMPING = .false.
@@ -190,74 +190,17 @@
end subroutine fault_setup
-
!==============================================================================================================
-! creates global indexing array ibool
-subroutine crm_ext_setup_indexing_fault(ibools, &
- xtemp,ytemp,ztemp,nspec,nglob,npointot)
- use generate_databases_par, only: nodes_coords_ext_mesh
+subroutine loading_coords_and_setup(faults,nglob,nspec,prname,myrank)
-! number of spectral elements in each block
- integer, intent(in) :: nspec,npointot
- integer, intent(out) :: nglob
-
-! arrays with the mesh
- integer, dimension(NGLLX,NGLLY,NGLLZ,nspec), intent(out) :: ibools
- double precision, dimension(NGLLX,NGLLY,NGLLZ,nspec), intent(in) :: xtemp,ytemp,ztemp
-
-! local parameters
-! variables for creating array ibools
- double precision, dimension(npointot) :: xp,yp,zp
- integer, dimension(npointot) :: locval
- logical, dimension(npointot) :: ifseg
-
- integer :: ieoff,ilocnum
- integer :: i,j,k,ispec
-
-! reshapes the arrays of GLL nodal coordinates into vectors
- do ispec=1,nspec
- ieoff = NGLLX * NGLLY * NGLLZ * (ispec-1)
- ilocnum = 0
- do k=1,NGLLZ
- do j=1,NGLLY
- do i=1,NGLLX
- ilocnum = ilocnum + 1
- xp(ilocnum+ieoff) = xtemp(i,j,k,ispec)
- yp(ilocnum+ieoff) = ytemp(i,j,k,ispec)
- zp(ilocnum+ieoff) = ztemp(i,j,k,ispec)
- enddo
- enddo
- enddo
- enddo
-
-! gets ibool indexing from local (GLL points) to global points
- call get_global(nspec,xp,yp,zp,ibools,locval,ifseg,nglob,npointot, &
- minval(nodes_coords_ext_mesh(1,:)),maxval(nodes_coords_ext_mesh(1,:)))
-
-!to do: try if the following works
-! call get_global(nspec,xtemp,ytemp,ztemp, ...
-! Fortran should automatically reshape xtemp into a vector. That's how ibools is passed.
-! If it works we don't need xp,yp,zp anymore.
-
-!- we can create a new indirect addressing to reduce cache misses
- call get_global_indirect_addressing(nspec,nglob,ibools)
-
-end subroutine crm_ext_setup_indexing_fault
-
-
-!==============================================================================================================
-
-subroutine loading_coords_and_setup(fdb,nglob,nspec,prname,myrank)
-
- type(fault_db_type), intent(inout) :: fdb
- integer, intent(in) :: nspec,nglob
+ integer, intent(in) :: nspec,nglob,i
integer, dimension(NGLLX,NGLLY,NGLLZ,nspec), intent(in) :: ibool
integer,intent(in) :: myrank
character(len=256), intent(in) :: prname ! 'proc***'
- integer :: nspec_fault,ifault
+ integer :: ifault,iflt
integer, dimension(3,NGLLSQUARE,nspec*6) :: ijk1, ijk2
integer :: ijk_face(3,NGLLX,NGLLY)
! The tolerance in get_global is SMALLVALTOTAL=1-10*(whole size of the model).
@@ -270,53 +213,45 @@
stop
endif
- read(IIN,*) nspec_fault
- allocate(fdb%ispec1(nspec_fault))
- allocate(fdb%iface1(nspec_fault))
- allocate(fdb%ijk1(3,NGLLX*NGLLY,nspec_fault))
+ do iflt=1,size(faults)
- allocate(fdb%ispec2(nspec_fault))
- allocate(fdb%iface2(nspec_fault))
- allocate(fdb%ijk2(3,NGLLX*NGLLY,nspec_fault))
+ read(IIN,*) faults(iflt)%nspec
+ if (faults(iflt)%nspec == 0) cycle
+ allocate(faults(iflt)%ispec1(fdb%nspec))
+ allocate(faults(iflt)%iface1(fdb%nspec))
+ allocate(faults(iflt)%ijk1(3,NGLLX*NGLLY,fdb%nspec))
+ allocate(faults(iflt)%ispec2(fdb%nspec))
+ allocate(faults(iflt)%iface2(fdb%nspec))
+ allocate(faults(iflt)%ijk2(3,NGLLX*NGLLY,fdb%nspec))
+ ! loading ispec1 ispec2 iface1 iface2 of fault elements.
+ do i=1,faults(iflt)%nspec
+ read(IIN,*) faults(iflt)%ispec1(i),faults(iflt)%ispec2(i),faults(iflt)%iface1(i),faults(iflt)%iface2(i)
+ enddo
- do i=1,nspec_fault
- read(IIN,*) fdb%ispec1(i),fdb%ispec2(i),fdb%iface1(i),fdb%iface2(i)
+ do ifault=1,faults(iflt)%nspec
+ ! we have identified a new fault element on fault side 1
+ iface_ref1 = faults(iflt)%iface1(ifault)
+ iface_ref2 = faults(iflt)%iface2(ifault)
+ ! gets i,j,k indices of GLL nodes in element face
+ call get_element_face_gll_indices(iface_ref1,ijk_face1,NGLLX,NGLLY)
+ call get_element_face_gll_indices(iface_ref2,ijk_face2,NGLLX,NGLLY)
+ igll = 0
+ do j=1,NGLLY
+ do i=1,NGLLX
+ igll = igll + 1
+ ijk1(:,igll,ifault)=ijk_face1(:,i,j) ! saving gll points of side 1 , needed for iulk1.
+ ijk2(:,igll,ifault)=ijk_face2(:,i,j) ! saving gll points of side 2 , needed for iulk1.
+ enddo
+ enddo
+ enddo
+ faults(iflt)%ijk1 = ijk1(:,:,1:faults(iflt)%nspec)
+ faults(iflt)%ijk2 = ijk2(:,:,1:faults(iflt)%nspec)
+
enddo
-
+
close(IIN)
- do ifault=1,nspec_fault
-
- ! we have identified a new fault element on fault side 1
- iface_ref1 = fdb%iface1(i)
- iface_ref2 = fdb%iface2(i)
-
- ! gets i,j,k indices of GLL nodes in element face
- call get_element_face_gll_indices(iface_ref1,ijk_face1,NGLLX,NGLLY)
- call get_element_face_gll_indices(iface_ref2,ijk_face2,NGLLX,NGLLY)
-
- igll = 0
- do j=1,NGLLY
- do i=1,NGLLX
-
- igll = igll + 1
-
- ijk1(:,igll,ifault)=ijk_face1(:,i,j) ! saving gll points of side 1 , needed for iulk1.
- ijk2(:,igll,ifault)=ijk_face2(:,i,j) ! saving gll points of side 2 , needed for iulk1.
-
- enddo
- enddo
- enddo
-
- fdb%ispec1 = ispec1(1:nspec_fault)
- fdb%iface1 = iface1(1:nspec_fault
- fdb%ijk1 = ijk1(:,:,1:nspec_fault)
-
- fdb%ispec2 = ispec2(1:nspec_fault)
- fdb%iface2 = iface2(1:nspec_fault)
- fdb%ijk2 = ijk2(:,:,1:nspec_fault)
-
end subroutine loading_coords_and_setup
!=============================================================================================================
More information about the CIG-COMMITS
mailing list