[cig-commits] r18663 - seismo/3D/FAULT_SOURCE/trunk/EXAMPLES/splay_faults
percygalvez at geodynamics.org
percygalvez at geodynamics.org
Tue Jun 28 05:00:53 PDT 2011
Author: percygalvez
Date: 2011-06-28 05:00:53 -0700 (Tue, 28 Jun 2011)
New Revision: 18663
Added:
seismo/3D/FAULT_SOURCE/trunk/EXAMPLES/splay_faults/save_fault_nodes_elements.py
Modified:
seismo/3D/FAULT_SOURCE/trunk/EXAMPLES/splay_faults/save_fault_crack_nodes.py
seismo/3D/FAULT_SOURCE/trunk/EXAMPLES/splay_faults/tsunami.py
Log:
adding save_fault_nodes_elements to save fault nodes (fault_down) reference
Modified: seismo/3D/FAULT_SOURCE/trunk/EXAMPLES/splay_faults/save_fault_crack_nodes.py
===================================================================
--- seismo/3D/FAULT_SOURCE/trunk/EXAMPLES/splay_faults/save_fault_crack_nodes.py 2011-06-24 14:25:25 UTC (rev 18662)
+++ seismo/3D/FAULT_SOURCE/trunk/EXAMPLES/splay_faults/save_fault_crack_nodes.py 2011-06-28 12:00:53 UTC (rev 18663)
@@ -1,54 +1,74 @@
########### Fault elements and nodes ###############
def (fault_id,fault_up,fault_down)
+ import cubit
+ import sys
+ import os
+ import numarray
- os.system('mkdir -p MESH')
+
+ os.system('mkdir -p MESH')
- # fault_u = 2 # fault surface up : surface 2
- # fault_d = 3 # fault surface down : surface 3
- txt =''
+ # fault_u # fault surface up
+ # fault_d # fault surface down
+
+ txt =''
+ filename = 'MESH/fault_file_'+str(fault_id)+'.dat'
- fault_file=open('MESH/fault_file_1.dat','w')
- list_hex=cubit.parse_cubit_list('hex','all')
-
- quads_fault_u = cubit.get_surface_quads(fault_u)
- quads_fault_d = cubit.get_surface_quads(fault_d)
-
- # TO DO : stop python properly in case fault nodes at both sides
- # do not match.
- if len(quads_fault_u) != len(quads_fault_d):
- stop
- # Writting number of elements at both sides to make
- # double sure everything is going fine .
- txt='%10i %10i\n' % (len(quads_fault_u),len(quads_fault_d))
- fault_file.write(txt)
-
- dic_quads_fault_u = dict(zip(quads_fault_u,quads_fault_u))
- dic_quads_fault_d = dict(zip(quads_fault_d,quads_fault_d))
-
- # FAULT SIDE UP
- for h in list_hex:
- faces = cubit.get_sub_elements('hex',h,2)
- for f in faces:
- if dic_quads_fault_u.has_key(f):
- nodes=cubit.get_connectivity('Face',f)
- # print 'h,fault nodes side up :',h,nodes[0],nodes[1],nodes[2],nodes[3]
- txt='%10i %10i %10i %10i %10i\n' % (h,nodes[0],\
- nodes[1],nodes[2],nodes[3])
- fault_file.write(txt)
-
- # FAULT SIDE DOWN
- for h in list_hex:
- faces = cubit.get_sub_elements('hex',h,2)
- for f in faces:
- if dic_quads_fault_d.has_key(f):
- nodes=cubit.get_connectivity('Face',f)
- # print 'h,fault nodes side down :',h,nodes[0],nodes[1],nodes[2],nodes[3]
- txt='%10i %10i %10i %10i %10i\n' % (h,nodes[0],\
- nodes[1],nodes[2],nodes[3])
- fault_file.write(txt)
-
- # CLOSING FAULT FILE
- fault_file.close()
-
+ fault_file=open(filename,'w')
+ list_hex=cubit.parse_cubit_list('hex','all')
+ quads_fault_u =[]
+ quads_fault_d =[]
+
+## Collecting upper fault surfaces
+ for fault in fault_up
+ quads_fault = cubit.get_surface_quads(fault)
+ quads_fault_u.append(quads_fault)
+
+## Collecting lower fault surfaces.
+ for fault in fault_down
+ quads_fault = cubit.get_surface_quads(fault)
+ quads_fault_d.append(quads_fault)
+
+ # TO DO : stop python properly in case fault nodes at both sides
+ # do not match.
+ if len(quads_fault_u) != len(quads_fault_d):
+ stop
+ # Writting number of elements at both sides to make
+ # double sure everything is going fine .
+ txt='%10i %10i\n' % (len(quads_fault_u),len(quads_fault_d))
+ fault_file.write(txt)
+
+ dic_quads_fault_u = dict(zip(quads_fault_u,quads_fault_u))
+ dic_quads_fault_d = dict(zip(quads_fault_d,quads_fault_d))
+
+
+## Looking at fault elements.
+
+ # FAULT SIDE DOWN
+ for h in list_hex:
+ faces = cubit.get_sub_elements('hex',h,2)
+ for f in faces:
+ if dic_quads_fault_d.has_key(f):
+ nodes=cubit.get_connectivity('Face',f)
+ # print 'h,fault nodes side down :',h,nodes[0],nodes[1],nodes[2],nodes[3]
+ txt='%10i %10i %10i %10i %10i\n' % (h,nodes[0],\
+ nodes[1],nodes[2],nodes[3])
+ fault_file.write(txt)
+
+ # FAULT SIDE UP
+ for h in list_hex:
+ faces = cubit.get_sub_elements('hex',h,2)
+ for f in faces:
+ if dic_quads_fault_u.has_key(f):
+ nodes=cubit.get_connectivity('Face',f)
+ # print 'h,fault nodes side up :',h,nodes[0],nodes[1],nodes[2],nodes[3]
+ txt='%10i %10i %10i %10i %10i\n' % (h,nodes[0],\
+ nodes[1],nodes[2],nodes[3])
+ fault_file.write(txt)
+
+ # CLOSING FAULT FILE
+ fault_file.close()
+
+
Added: seismo/3D/FAULT_SOURCE/trunk/EXAMPLES/splay_faults/save_fault_nodes_elements.py
===================================================================
--- seismo/3D/FAULT_SOURCE/trunk/EXAMPLES/splay_faults/save_fault_nodes_elements.py (rev 0)
+++ seismo/3D/FAULT_SOURCE/trunk/EXAMPLES/splay_faults/save_fault_nodes_elements.py 2011-06-28 12:00:53 UTC (rev 18663)
@@ -0,0 +1,69 @@
+#!python
+# Opening fault cracks
+# Input : surface up and down.
+import cubit
+
+class fault_input:
+ def __init__(self,id,surface_u,surface_d):
+ self.id = id
+ self.surface_u = surface_u
+ self.surface_d = surface_d
+ self.name = 'MESH/fault_file_'+str(id)+'.dat'
+
+def save_cracks(name,list_surface_up,list_surface_down):
+ quads_fault_up = []
+ quads_fault_down = []
+ for surface in list_surface_up :
+ quads_fault = cubit.get_surface_quads(surface)
+ quads_fault_up.append(quads_fault)
+ for surface in list_surface_down :
+ quads_fault = cubit.get_surface_quads(surface)
+ quads_fault_down.append(quads_fault)
+ # TO DO : stop python properly in case fault nodes at both sides
+ # do not match.
+ # if len(quads_fault_u) != len(quads_fault_d): stop
+ #
+ # SAVING FAULT ELEMENTS AND NODES
+ return quads_fault_up,quads_fault_down
+
+def unpack_list(fault_list):
+ list_fault = []
+ for i in range(0,len(fault_list)):
+ el=list(fault_list[i])
+ for j in el:
+ list_fault.append(j)
+ return list_fault
+
+def save_elements_nodes(name,quads_fault_u,quads_fault_d):
+ fault_file = open(name,'w')
+ txt =''
+ list_hex=cubit.parse_cubit_list('hex','all')
+ txt='%10i %10i\n' % (len(quads_fault_u),len(quads_fault_d))
+ fault_file.write(txt)
+
+ dic_quads_fault_u = dict(zip(quads_fault_u,quads_fault_u))
+ dic_quads_fault_d = dict(zip(quads_fault_d,quads_fault_d))
+
+ # FAULT SIDE DOWN
+ for h in list_hex:
+ faces = cubit.get_sub_elements('hex',h,2)
+ for f in faces:
+ if dic_quads_fault_d.has_key(f):
+ nodes=cubit.get_connectivity('Face',f)
+ # print 'h,fault nodes side down :',h,nodes[0],nodes[1],nodes[2],nodes[3]
+ txt='%10i %10i %10i %10i %10i\n' % (h,nodes[0],\
+ nodes[1],nodes[2],nodes[3])
+ fault_file.write(txt)
+
+ # FAULT SIDE UP
+ for h in list_hex:
+ faces = cubit.get_sub_elements('hex',h,2)
+ for f in faces:
+ if dic_quads_fault_u.has_key(f):
+ nodes=cubit.get_connectivity('Face',f)
+ # print 'h,fault nodes side up :',h,nodes[0],nodes[1],nodes[2],nodes[3]
+ txt='%10i %10i %10i %10i %10i\n' % (h,nodes[0],\
+ nodes[1],nodes[2],nodes[3])
+ fault_file.write(txt)
+
+ fault_file.close()
Modified: seismo/3D/FAULT_SOURCE/trunk/EXAMPLES/splay_faults/tsunami.py
===================================================================
--- seismo/3D/FAULT_SOURCE/trunk/EXAMPLES/splay_faults/tsunami.py 2011-06-24 14:25:25 UTC (rev 18662)
+++ seismo/3D/FAULT_SOURCE/trunk/EXAMPLES/splay_faults/tsunami.py 2011-06-28 12:00:53 UTC (rev 18663)
@@ -12,8 +12,8 @@
km = 1000
y_vert=0*km
-h = 1
-h1 = 0.1
+h = 0 #0.1#(h=1)
+h1 = 0 #0.1#(h1=0.1)
# dipping angle 1
d = math.pi/180
dip1=7.5*d
@@ -35,8 +35,8 @@
# Bulk points ###
xbulk.append(-10*km) #x1
-xbulk.append(250*km) #x2
-xbulk.append(250*km) #x3
+xbulk.append(240*km) #x2
+xbulk.append(240*km) #x3
xbulk.append(-10*km) #x4
zbulk.append(0*km) #y1
@@ -48,42 +48,28 @@
### CRACKS ##########
x.append(0*km) #x5
-x.append(4*km/math.tan(dip1)) #x6
+x.append(4*km/math.tan(dip1)) #x6 = x11 (triple joint)
x.append(x[1]+8*km/math.tan(dip2)) #x7
-x.append(90*km) #x8
-x.append(227*km) #x9
-x.append(x[1]) #x10
+x.append(x[2]) #x8 = x7 (kink point)
+x.append(90*km) #x9
+x.append(90*km+h1) #x10 = x9+h1 (surface point, fault C)
+x.append(x[1]) #x11 = x6
+x.append(227*km) #x12
+x.append(x[7]+h1) #x13 = x12+h1 (surface point, fault B)
-z.append(-30*km) #z5
-z.append(-26*km+h1) #z6 ! shifting dh=h1 joint points.
-z.append(-18*km) #z7
-z.append(0) #z8
-z.append(0) #z9
-z.append(-26*km-h1) #z10
+z.append(-30*km) #z5
+z.append(-26*km) #z6
+z.append(-18*km) #z7 (kink point)
+z.append(z[2]-h1) #z8 = z7-h1 (kink point)
+z.append(0) #z9
+z.append(0) #z10 = z9
+z.append(-26*km) #z11 = z6
+z.append(0) #z12
+z.append(0) #z13
-y=[y_vert]*6
+y=[y_vert]*9
-### 6 midpoints for three fault cracks A, B and C ###
-# Fault A (already open)
-#xmid.append((x[1])) #9 (52) ! JOINT POINT
-#zmid.append((z[1]))
-#xmid.append((x[5])) #10 (52) ! JOINT POINT
-#zmid.append((z[5]))
-
-# Fault B
-xmid.append((x[2]+x[1])/2) #11 (76)
-xmid.append((x[2]+x[1])/2) #12 (76)
-zmid.append((z[2]+z[1]+h)/2)
-zmid.append((z[2]+z[1]-h)/2)
-
-# Fault C
-xmid.append((x[3]+x[2])/2) #13 (87)
-xmid.append((x[3]+x[2])/2) #14 (87)
-zmid.append((z[3]+z[2]+h)/2)
-zmid.append((z[3]+z[2]-h)/2)
-
-
ymid=[y_vert]*6
#####
@@ -106,12 +92,12 @@
################ creating fault domains #################################
-bulk1="create curve vertex 1 8" #c1
-bulk2="create curve vertex 8 9" #c2
-bulk3="create curve vertex 9 2" #c3
-bulk4="create curve vertex 2 3" #c4
-bulk5="create curve vertex 3 4" #c5
-bulk6="create curve vertex 4 1" #c6
+bulk1="create curve vertex 1 9" #c1
+bulk2="create curve vertex 10 12" #c2
+bulk3="create curve vertex 13 2" #c3
+bulk4="create curve vertex 2 3" #c4
+bulk5="create curve vertex 3 4" #c5
+bulk6="create curve vertex 4 1" #c6
cubit.cmd(bulk1)
cubit.cmd(bulk2)
@@ -120,63 +106,96 @@
cubit.cmd(bulk5)
cubit.cmd(bulk6)
-fault_up_A1="create curve spline vertex 5 6" #c7
-fault_up_A2="create curve spline vertex 6 9" #c8
+fault_up_A1="create curve spline vertex 5 6" #c7
+fault_up_A2="create curve spline vertex 6 12" #c8
-fault_down_A1="create curve spline vertex 5 10" #c9
-fault_down_A2="create curve spline vertex 10 9" #c10
+fault_down_A1="create curve spline vertex 5 11" #c9
+fault_down_A2="create curve spline vertex 11 13" #c10
-fault_up_B="create curve spline vertex 6 11 7" #c11
-fault_down_B="create curve spline vertex 6 12 7" #c12
+fault_up_BC1="create curve vertex 9 7" #c11
+fault_up_BC2="create curve vertex 7 6" #c12
+fault_down_BC1="create curve vertex 10 8" #c13
+fault_down_BC2="create curve vertex 8 6" #c14
-fault_up_C="create curve spline vertex 7 13 8" #c13
-fault_down_C="create curve spline vertex 7 14 8" #c14
-
cubit.cmd(fault_up_A1)
cubit.cmd(fault_up_A2)
cubit.cmd(fault_down_A1)
cubit.cmd(fault_down_A2)
-cubit.cmd(fault_up_B)
-cubit.cmd(fault_down_B)
-cubit.cmd(fault_up_C)
-cubit.cmd(fault_down_C)
-surface="create surface curve 1 13 11 7 9 10 3 4 5 6 "
+cubit.cmd(fault_up_BC1)
+cubit.cmd(fault_up_BC2)
+cubit.cmd(fault_down_BC1)
+cubit.cmd(fault_down_BC2)
+
+#surface="create surface curve 1 11 12 7 9 10 3 4 5 6"
+surface="create surface curve 1 11 12 10 3 4 5 6"
+
cubit.cmd(surface)
-surface="create surface curve 2 14 12 8"
+surface="create surface curve 2 13 14 8"
cubit.cmd(surface)
+
cubit.cmd("sweep surface 1 vector 0 1 0 distance "+str(90*km))
-cubit.cmd("sweep surface 2 vector 0 1 0 distance "+str(90*km))
+cubit.cmd("sweep surface 2 vector 0 1 0 distance "+str(90*km))
+cubit.cmd("sweep curve 7 vector 0 1 0 distance "+str(90*km))
+cubit.cmd("sweep curve 9 vector 0 1 0 distance "+str(90*km))
+# fault_up_A1
+cubit.cmd("sweep curve 7 vector 0 1 0 distance")
+#fault_down_A1
+cubit.cmd("sweep curve 9 vector 0 1 0 distance")
+### Coarsening ########
+xcoorse=[]
+ycoorse=[]
+zcoorse=[]
-### fault crack (Not necessary here) ###
-# FAULT A
-#cubit.cmd('curve 8 7 9 10 merge off')
-# FAULT B
-#cubit.cmd('curve 11 12 merge off')
-# FAULT C
-#cubit.cmd('curve 13 14 merge off')
+kx=250*km
+ky=90*km
+kz=50*km
-#cubit.cmd('merge tolerance 1e-3')
-###
+# FAULT_A1_up
+#cubit.cmd("curve 7 merge off")
+#cubit.cmd("curve 24 merge off")
+cubit.cmd("surface 3 merge off ")
+# FAULT_A1_down
+#cubit.cmd("curve 9 merge off")
+#cubit.cmd("curve 26 merge off")
+cubit.cmd("surface 3 merge off ")
-#####################################################
-elementsize = 2500
-
# IMPRINTING
cubit.cmd("imprint all")
-
# MERGING
cubit.cmd("merge all")
+# Sticking fault borders ##
+#### Fault A.
+# UP
+cubit.cmd("merge curve 8 with curve 10 force") # init
+cubit.cmd("merge curve 22 with curve 34 force") # final
+## DOWN.
+#cubit.cmd("merge curve 7 with curve 9 force")
+#cubit.cmd('merge curve 24 with curve 26 force')
+##### Fault BC
+#### Up
+cubit.cmd("merge curve 11 with curve 13 force")
+cubit.cmd("merge curve 30 with curve 42 force")
+#### DOWN
+cubit.cmd("merge curve 12 with curve 14 force")
+cubit.cmd("merge curve 28 with curve 35 force")
-cubit.cmd("surface 13 18 size "+str(elementsize))
-cubit.cmd("volume 1 2 3 size "+str(elementsize))
-cubit.cmd("surface 13 18 scheme pave")
-cubit.cmd("mesh surface 13 18 ")
-cubit.cmd("mesh volume 1 2 3")
+# Mesh coinciding fault_A upper and lower boundaries .
+# Meshing the whole domain.
+#cubit.cmd("surface 13 18 size "+str(elementsize))
+#cubit.cmd("volume 1 2 3 size "+str(elementsize))
+#cubit.cmd("surface 13 18 scheme pave")
+#cubit.cmd("mesh surface 13 18 ")
+#cubit.cmd("mesh volume 1 2 3")
#cubit.cmd("unmerge surface 2 3")
-########## Cracks #######
+########## Coarse shield #######
+
+
+
+
+
More information about the CIG-COMMITS
mailing list