[cig-commits] r18335 - in seismo/3D/FAULT_SOURCE/branches/new_fault_db: CUBIT decompose_mesh_SCOTCH src
percygalvez at geodynamics.org
percygalvez at geodynamics.org
Mon May 9 08:41:36 PDT 2011
Author: percygalvez
Date: 2011-05-09 08:41:36 -0700 (Mon, 09 May 2011)
New Revision: 18335
Added:
seismo/3D/FAULT_SOURCE/branches/new_fault_db/CUBIT/boundary_definition.py
seismo/3D/FAULT_SOURCE/branches/new_fault_db/CUBIT/cubit2specfem3d.py
seismo/3D/FAULT_SOURCE/branches/new_fault_db/CUBIT/fault_elements_1.dat
Modified:
seismo/3D/FAULT_SOURCE/branches/new_fault_db/CUBIT/crack_fault.py
seismo/3D/FAULT_SOURCE/branches/new_fault_db/decompose_mesh_SCOTCH/fault_scotch.f90
seismo/3D/FAULT_SOURCE/branches/new_fault_db/src/fault_object.f90
Log:
PYTHON/CUBIT fault_elements directly
Added: seismo/3D/FAULT_SOURCE/branches/new_fault_db/CUBIT/boundary_definition.py
===================================================================
--- seismo/3D/FAULT_SOURCE/branches/new_fault_db/CUBIT/boundary_definition.py (rev 0)
+++ seismo/3D/FAULT_SOURCE/branches/new_fault_db/CUBIT/boundary_definition.py 2011-05-09 15:41:36 UTC (rev 18335)
@@ -0,0 +1,394 @@
+#############################################################################
+# boundary_definition.py #
+# this file is part of GEOCUBIT #
+# #
+# Created by Emanuele Casarotti #
+# Copyright (c) 2008 Istituto Nazionale di Geofisica e Vulcanologia #
+# #
+#############################################################################
+# #
+# GEOCUBIT is free software: you can redistribute it and/or modify #
+# it under the terms of the GNU General Public License as published by #
+# the Free Software Foundation, either version 3 of the License, or #
+# (at your option) any later version. #
+# #
+# GEOCUBIT is distributed in the hope that it will be useful, #
+# but WITHOUT ANY WARRANTY; without even the implied warranty of #
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the #
+# GNU General Public License for more details. #
+# #
+# You should have received a copy of the GNU General Public License #
+# along with GEOCUBIT. If not, see <http://www.gnu.org/licenses/>. #
+# #
+#############################################################################
+
+def define_absorbing_surf():
+ """
+ define the absorbing surfaces for a layered topological box where boundary are surfaces parallel to the axis.
+ it returns absorbing_surf,absorbing_surf_xmin,absorbing_surf_xmax,absorbing_surf_ymin,absorbing_surf_ymax,absorbing_surf_bottom,topo_surf
+ where
+ absorbing_surf is the list of all the absorbing boundary surf
+ absorbing_surf_xmin is the list of the absorbing boundary surfaces that correnspond to x=xmin
+ ...
+ absorbing_surf_bottom is the list of the absorbing boundary surfaces that correspond to z=zmin
+ """
+ try:
+ cubit.cmd('comment')
+ except:
+ try:
+ import cubit
+ cubit.init([""])
+ except:
+ print 'error importing cubit'
+ import sys
+ sys.exit()
+ absorbing_surf=[]
+ absorbing_surf_xmin=[]
+ absorbing_surf_xmax=[]
+ absorbing_surf_ymin=[]
+ absorbing_surf_ymax=[]
+ absorbing_surf_bottom=[]
+ top_surf=[]
+
+
+ list_vol=cubit.parse_cubit_list("volume","all")
+ init_n_vol=len(list_vol)
+ zmax_box=cubit.get_total_bounding_box("volume",list_vol)[7]
+ zmin_box=cubit.get_total_bounding_box("volume",list_vol)[6] #it is the z_min of the box ... box= xmin,xmax,d,ymin,ymax,d,zmin...
+ xmin_box=cubit.get_total_bounding_box("volume",list_vol)[0]
+ xmax_box=cubit.get_total_bounding_box("volume",list_vol)[1]
+ ymin_box=cubit.get_total_bounding_box("volume",list_vol)[3]
+ ymax_box=cubit.get_total_bounding_box("volume",list_vol)[4]
+ list_surf=cubit.parse_cubit_list("surface","all")
+
+ print '##boundary box: '
+ print '## x min: ' + str(xmin_box)
+ print '## y min: ' + str(ymin_box)
+ print '## z min: ' + str(zmin_box)
+ print '## x max: ' + str(xmax_box)
+ print '## y max: ' + str(ymax_box)
+ print '## z max: ' + str(zmax_box)
+
+
+
+
+# for k in list_surf:
+# center_point = cubit.get_center_point("surface", k)
+# if abs((center_point[0] - xmin_box)/xmin_box) <= 0.005:
+# absorbing_surf_xmin.append(k)
+# absorbing_surf.append(k)
+# elif abs((center_point[0] - xmax_box)/xmax_box) <= 0.005:
+# absorbing_surf_xmax.append(k)
+# absorbing_surf.append(k)
+# elif abs((center_point[1] - ymin_box)/ymin_box) <= 0.005:
+# absorbing_surf_ymin.append(k)
+# absorbing_surf.append(k)
+# elif abs((center_point[1] - ymax_box)/ymax_box) <= 0.005:
+# absorbing_surf_ymax.append(k)
+# absorbing_surf.append(k)
+# elif abs((center_point[2] - zmin_box)/zmin_box) <= 0.005:
+# absorbing_surf_bottom.append(k)
+# absorbing_surf.append(k)
+# else:
+# sbox=cubit.get_bounding_box('surface',k)
+# dz=abs((sbox[7] - zmax_box)/zmax_box)
+# normal=cubit.get_surface_normal(k)
+# zn=normal[2]
+# dn=abs(zn-1)
+# if dz <= 0.001 and dn < 0.2:
+# top_surf.append(k)
+
+ #box lengths
+ x_len = abs( xmax_box - xmin_box)
+ y_len = abs( ymax_box - ymin_box)
+ z_len = abs( zmax_box - zmin_box)
+
+ print '##boundary box: '
+ print '## x length: ' + str(x_len)
+ print '## y length: ' + str(y_len)
+ print '## z length: ' + str(z_len)
+
+ # tolerance parameters
+ absorbing_surface_distance_tolerance=0.005
+ topographic_surface_distance_tolerance=0.001
+ topographic_surface_normal_tolerance=0.2
+
+ for k in list_surf:
+ center_point = cubit.get_center_point("surface", k)
+ if abs((center_point[0] - xmin_box)/x_len) <= absorbing_surface_distance_tolerance:
+ absorbing_surf_xmin.append(k)
+ absorbing_surf.append(k)
+ elif abs((center_point[0] - xmax_box)/x_len) <= absorbing_surface_distance_tolerance:
+ absorbing_surf_xmax.append(k)
+ absorbing_surf.append(k)
+ elif abs((center_point[1] - ymin_box)/y_len) <= absorbing_surface_distance_tolerance:
+ absorbing_surf_ymin.append(k)
+ absorbing_surf.append(k)
+ elif abs((center_point[1] - ymax_box)/y_len) <= absorbing_surface_distance_tolerance:
+ absorbing_surf_ymax.append(k)
+ absorbing_surf.append(k)
+ elif abs((center_point[2] - zmin_box)/z_len) <= absorbing_surface_distance_tolerance:
+ print 'center_point[2]' + str(center_point[2])
+ print 'kz:' + str(k)
+ absorbing_surf_bottom.append(k)
+ absorbing_surf.append(k)
+
+
+ else:
+ sbox=cubit.get_bounding_box('surface',k)
+ dz=abs((sbox[7] - zmax_box)/z_len)
+ normal=cubit.get_surface_normal(k)
+ zn=normal[2]
+ dn=abs(zn-1)
+ if dz <= topographic_surface_distance_tolerance and dn < topographic_surface_normal_tolerance:
+ top_surf.append(k)
+
+ return absorbing_surf,absorbing_surf_xmin,absorbing_surf_xmax,absorbing_surf_ymin,absorbing_surf_ymax,absorbing_surf_bottom,top_surf
+
+def define_absorbing_surf_nopar():
+ """
+ define the absorbing surfaces for a layered topological box where boundary surfaces are not parallel to the axis.
+ it returns absorbing_surf,topo_surf
+ where
+ absorbing_surf is the list of all the absorbing boundary surf
+ """
+ try:
+ cubit.cmd('comment')
+ except:
+ try:
+ import cubit
+ cubit.init([""])
+ except:
+ print 'error importing cubit'
+ import sys
+ sys.exit()
+ from sets import Set
+ def product(*args, **kwds):
+ # product('ABCD', 'xy') --> Ax Ay Bx By Cx Cy Dx Dy
+ # product(range(2), repeat=3) --> 000 001 010 011 100 101 110 111
+ pools = map(tuple, args) * kwds.get('repeat', 1)
+ result = [[]]
+ for pool in pools:
+ result = [x+[y] for x in result for y in pool]
+ return result
+ absorbing_surf=[]
+ absorbing_surf_xmin=[]
+ absorbing_surf_xmax=[]
+ absorbing_surf_ymin=[]
+ absorbing_surf_ymax=[]
+ absorbing_surf_bottom=[]
+ top_surf=[]
+ list_vol=cubit.parse_cubit_list("volume","all")
+ init_n_vol=len(list_vol)
+ zmax_box=cubit.get_total_bounding_box("volume",list_vol)[7]
+ zmin_box=cubit.get_total_bounding_box("volume",list_vol)[6] #it is the z_min of the box ... box= xmin,xmax,d,ymin,ymax,d,zmin...
+ xmin_box=cubit.get_total_bounding_box("volume",list_vol)[0]
+ xmax_box=cubit.get_total_bounding_box("volume",list_vol)[1]
+ ymin_box=cubit.get_total_bounding_box("volume",list_vol)[3]
+ ymax_box=cubit.get_total_bounding_box("volume",list_vol)[4]
+ list_surf=cubit.parse_cubit_list("surface","all")
+ lv=[]
+ for k in list_surf:
+ sbox=cubit.get_bounding_box('surface',k)
+ dzmax=abs((sbox[7] - zmax_box)/zmax_box)
+ dzmin=abs((sbox[6] - zmin_box)/zmin_box)
+ normal=cubit.get_surface_normal(k)
+ zn=normal[2]
+ if dzmax <= 0.001 and zn > 0.7:
+ top_surf.append(k)
+ list_vertex=cubit.get_relatives('surface',k,'vertex')
+ for v in list_vertex:
+ valence=cubit.get_valence(v)
+ if valence <= 4: #valence 3 is a corner, 4 is a vertex between 2 volumes, > 4 is a vertex not in the boundaries
+ lv.append(v)
+ elif dzmin <= 0.001 and zn < -0.7:
+ absorbing_surf.append(k)
+ lp=[]
+ combs=product(lv,lv)
+ for comb in combs:
+ v1=comb[0]
+ v2=comb[1]
+ c=Set(cubit.get_relatives("vertex",v1,"curve")) & Set(cubit.get_relatives("vertex",v2,"curve"))
+ if len(c) == 1:
+ p=cubit.get_center_point("curve",list(c)[0])
+ lp.append(p)
+ for k in list_surf:
+ center_point = cubit.get_center_point("surface", k)
+ for p in lp:
+ if abs((center_point[0] - p[0])/p[0]) <= 0.005 and abs((center_point[1] - p[1])/p[1]) <= 0.005:
+ absorbing_surf.append(k)
+ break
+ return absorbing_surf,top_surf
+
+def define_absorbing_surf_sphere():
+ try:
+ cubit.cmd('comment')
+ except:
+ try:
+ import cubit
+ cubit.init([""])
+ except:
+ print 'error importing cubit'
+ import sys
+ sys.exit()
+ surf=[]
+ list_surf=cubit.parse_cubit_list("surface","all")
+ for s in list_surf:
+ v=cubit.get_relatives('surface',s,'volume')
+ if len(v) == 1:
+ surf.append(s)
+ return surf
+
+def define_block():
+ try:
+ cubit.cmd('comment')
+ except:
+ try:
+ import cubit
+ cubit.init([""])
+ except:
+ print 'error importing cubit'
+ import sys
+ sys.exit()
+ list_vol=cubit.parse_cubit_list("volume","all")
+ init_n_vol=len(list_vol)
+ list_name=map(lambda x: 'vol'+x,map(str,list_vol))
+ return list_vol,list_name
+
+def build_block(vol_list,name):
+ from sets import Set
+ try:
+ cubit.cmd('comment')
+ except:
+ try:
+ import cubit
+ cubit.init([""])
+ except:
+ print 'error importing cubit'
+ import sys
+ sys.exit()
+ block_list=cubit.get_block_id_list()
+ if len(block_list) > 0:
+ id_block=max(block_list)
+ else:
+ id_block=0
+ for v,n in zip(vol_list,name):
+ id_block+=1
+ v_other=Set(vol_list)-Set([v])
+ #command= 'block '+str(id_block)+' hex in node in vol '+str(v)+' except hex in vol '+str(list(v_other))
+ command= 'block '+str(id_block)+' hex in vol '+str(v)
+ command = command.replace("["," ").replace("]"," ")
+ cubit.cmd(command)
+ command = "block "+str(id_block)+" name '"+n+"'"
+ cubit.cmd(command)
+
+def build_block_side(surf_list,name,obj='surface'):
+ try:
+ cubit.cmd('comment')
+ except:
+ try:
+ import cubit
+ cubit.init([""])
+ except:
+ print 'error importing cubit'
+ import sys
+ sys.exit()
+ id_nodeset=cubit.get_next_nodeset_id()
+ id_block=cubit.get_next_block_id()
+
+
+ if obj == 'hex':
+ txt='hex in node in surface'
+ txt1='block '+str(id_block)+ ' '+ txt +' '+str(list(surf_list))
+ txt2="block "+str(id_block)+" name '"+name+"'"
+ txt1=txt1.replace("["," ").replace("]"," ")
+ elif obj == 'node':
+ txt=obj+' in surface'
+ txt1= 'nodeset '+str(id_nodeset)+ ' '+ txt +' '+str(list(surf_list))
+ txt1 = txt1.replace("["," ").replace("]"," ")
+ txt2 = "nodeset "+str(id_nodeset)+" name '"+name+"'"
+ elif obj == 'face' or obj == 'edge':
+ txt=obj+' in surface'
+ txt1= 'block '+str(id_block)+ ' '+ txt +' '+str(list(surf_list))
+ txt1 = txt1.replace("["," ").replace("]"," ")
+ txt2 = "block "+str(id_block)+" name '"+name+"'"
+ else:
+ txt1=''
+ # do not execute: block id might be wrong
+ print "##block "+str(id_block)+" name '"+name+"_notsupported (only hex,face,edge,node)'"
+ txt2=''
+
+
+ cubit.cmd(txt1)
+ cubit.cmd(txt2)
+
+def define_bc(*args,**keys):
+ parallel=keys.get('parallel',True)
+ closed=keys.get('closed',False)
+ if not closed:
+ print "##open region"
+
+ if parallel:
+ surf,xmin,xmax,ymin,ymax,bottom,topo=define_absorbing_surf()
+ else:
+ surf,topo=define_absorbing_surf_nopar()
+ v_list,name_list=define_block()
+ build_block(v_list,name_list)
+ entities=args[0]
+ print entities
+ for entity in entities:
+ print "##entity: "+str(entity)
+
+ build_block_side(topo,entity+'_topo',obj=entity)
+ build_block_side(surf,entity+'_abs',obj=entity)
+ if parallel:
+ build_block_side(xmin,entity+'_abs_xmin',obj=entity)
+ build_block_side(xmax,entity+'_abs_xmax',obj=entity)
+ build_block_side(ymin,entity+'_abs_ymin',obj=entity)
+ build_block_side(ymax,entity+'_abs_ymax',obj=entity)
+ build_block_side(bottom,entity+'_abs_bottom',obj=entity)
+ else:
+ print "##closed region"
+
+ surf=define_absorbing_surf_sphere()
+ v_list,name_list=define_block()
+ build_block(v_list,name_list)
+ entities=args[0]
+ for entity in entities:
+ build_block_side(surf,entity+'_closedvol',obj=entity)
+
+
+
+#entities=['surface','face','edge','node','hex']
+#define_bc(entities,parallel=True)
+#define_bc(entities,parallel=False)
+#define_bc(entities,parallel=False,closed=True)
+
+
+# call
+#entities=['surface','face']
+#define_bc(entities,parallel=True)
+
+#block 1 attribute count 5
+#block 2 attribute count 0
+#block 2 name '/prova/interface1'
+#block 2 attribute count 3
+#block 3 attribute count 5
+#block 1 attribute index 2 1500
+#block 1 attribute count 2
+#block 3 attribute index 1 2
+#block 3 attribute index 2 5800
+#block 3 attribute index 3 3900
+#block 3 attribute index 4 1500
+#block 3 attribute index 5 3.5
+#block 1 name 'top'
+#block 3 name 'bottom'
+#block 2 attribute index 1 -1
+#block 2 attribute index 2 2
+#block 2 attribute count 4
+#block 2 attribute index 2 1
+#block 2 attribute index 3 2
+
+
+
+
Modified: seismo/3D/FAULT_SOURCE/branches/new_fault_db/CUBIT/crack_fault.py
===================================================================
--- seismo/3D/FAULT_SOURCE/branches/new_fault_db/CUBIT/crack_fault.py 2011-05-08 08:12:31 UTC (rev 18334)
+++ seismo/3D/FAULT_SOURCE/branches/new_fault_db/CUBIT/crack_fault.py 2011-05-09 15:41:36 UTC (rev 18335)
@@ -1,6 +1,8 @@
#!/usr/bin/env python
import cubit
+import boundary_definition
+import cubit2specfem3d
import math
import os
import sys
@@ -84,66 +86,86 @@
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)
+########### SIDESETS (NOT USED ) ###############
+#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)
-#### SIDESETS
+########### Fault elements and nodes ###############
-sidesetlist = cubit.get_sideset_id_list()
+os.system('mkdir -p MESH')
-for sidesetid in sidesetlist:
- sidesetname=cubit.get_exodus_entity_name('sideset',sidesetid)
- print 'sideset '+sidesetname+' ids'
- print 'sideset '+str(sidesetid)
+fault_u = 2 # fault surface up : surface 2
+fault_d = 3 # fault surface down : surface 3
+txt =''
+fault_file=open('MESH/fault_file_1.dat','w')
+list_hex=cubit.parse_cubit_list('hex','all')
-for sidesetid in sidesetlist:
- sidesetquads = cubit.get_sideset_quads(sidesetid)
- sidesetsurfaces= cubit.get_sideset_surfaces(sidesetid)
- print 'sideset quads :' , sidesetquads # empty no quads ?
- print 'sideset surface :' +str(sidesetsurfaces)
+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)
-# sidesetdescription = cubit.get_exodus_entity_description('sideset',sidesetid)
-# sidesetel = cubit.get_exodus_id('sideset',sidesetid)
-# face = cubit.get_sub_elements('sideset',sidesetid,2)
-# print 'sideset description :' + sidesetdescription # ....
-# print 'sideset el:' +str(sidesetel) # -1
-# print 'face : '+str(face)
+dic_quads_fault_u = dict(zip(quads_fault_u,quads_fault_u))
+dic_quads_fault_d = dict(zip(quads_fault_d,quads_fault_d))
-# sideset_list = cubit.parse_cubit_list('node',sidesetid)
-# name=cubit.get_exodus_entity_name("sideset",nset)
-# id_sideset=cubit.get_exodus_id("hex",nset)
-# print name
-# print id_sideset
+# 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 'nodes :',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)
-### INTERFACE BLOCK
-fault_u = 300
-fault_d = 301
-side_u = "block "+str(fault_u)+" surface 2"
-side_d = "block "+str(fault_d)+" surface 3"
-cubit.cmd(side_u)
-cubit.cmd(side_d)
-el_u=cubit.get_block_faces(fault_u)
-el_d=cubit.get_block_faces(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 'nodes :',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)
-for el in el_u:
- print 'el_u :', el
+# CLOSING FAULT FILE
+fault_file.close()
-for el in el_d:
- print 'el_d :', el
+# FOR THE BULK (Seismic wave propagation part for SESAME)
-#cubit.get_block_hexes(block_id)
+###### This is boundary_definition.py of GEOCUBIT
+#..... which extracts the bounding faces and defines them into blocks
+boundary_definition.entities=['face']
+boundary_definition.define_bc(boundary_definition.entities,parallel=True)
+
+#### Define material properties for the 2 volumes ################
+cubit.cmd('#### DEFINE MATERIAL PROPERTIES #######################')
+
+# Material properties in concordance with tpv5 benchmark.
+
+cubit.cmd('block 1 name "elastic 1" ') # material region
+cubit.cmd('block 1 attribute count 5')
+cubit.cmd('block 1 attribute index 1 1') # flag for fault side 1
+cubit.cmd('block 1 attribute index 2 6000') # vp
+cubit.cmd('block 1 attribute index 3 3464') # vs
+cubit.cmd('block 1 attribute index 4 2670') # rho
+cubit.cmd('block 1 attribute index 5 13') # Q flag (see constants.h: IATTENUATION_ ... )
+#### Export to SESAME format using cubit2specfem3d.py of GEOCUBIT
+
+cubit2specfem3d.export2SESAME('MESH')
+
+# all files needed by SCOTCH are now in directory MESH
-
-
-
-
-
-
-
-
Added: seismo/3D/FAULT_SOURCE/branches/new_fault_db/CUBIT/cubit2specfem3d.py
===================================================================
--- seismo/3D/FAULT_SOURCE/branches/new_fault_db/CUBIT/cubit2specfem3d.py (rev 0)
+++ seismo/3D/FAULT_SOURCE/branches/new_fault_db/CUBIT/cubit2specfem3d.py 2011-05-09 15:41:36 UTC (rev 18335)
@@ -0,0 +1,775 @@
+#!python
+#############################################################################
+# cubit2specfem3d.py #
+# this file is part of GEOCUBIT #
+# #
+# Created by Emanuele Casarotti #
+# Copyright (c) 2008 Istituto Nazionale di Geofisica e Vulcanologia #
+# #
+#############################################################################
+# #
+# GEOCUBIT is free software: you can redistribute it and/or modify #
+# it under the terms of the GNU General Public License as published by #
+# the Free Software Foundation, either version 3 of the License, or #
+# (at your option) any later version. #
+# #
+# GEOCUBIT is distributed in the hope that it will be useful, #
+# but WITHOUT ANY WARRANTY; without even the implied warranty of #
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the #
+# GNU General Public License for more details. #
+# #
+# You should have received a copy of the GNU General Public License #
+# along with GEOCUBIT. If not, see <http://www.gnu.org/licenses/>. #
+# #
+#############################################################################
+#
+#for a complete definition of the format of the mesh in SPECFEM3D_SESAME check the manual (http:/):
+#
+#USAGE
+#
+#############################################################################
+#PREREQUISITE
+#The mesh must be prepared
+# automatically using the module boundary_definition (see boundary_definition.py for more information)
+#or
+# manually following the convention:
+# - each material should have a block defined by:
+# material domain_flag (acoustic/elastic/poroelastic)name,flag of the material (integer),p velocity
+# (or the full description: name, flag, vp, vs, rho, Q ... if not present these last 3 parameters will be
+# interpolated by module mat_parameter)
+# - each mesh should have the block definition for the face on the free_surface (topography),
+# the name of this block must be 'face_topo' or you can change the default name in mesh.topo defined in profile.
+# - each mesh should have the block definition for the faces on the absorbing boundaries,
+# one block for each surface with x=Xmin,x=Xmax,y=Ymin,y=Ymax and z=bottom. The names of
+# the blocks should contain the strings "xmin,xmax,ymin,ymax,bottom"
+#
+#############################################################################
+#RUN
+#In a python script or in the cubit python tab call:
+#
+# export2SESAME(path_exporting_mesh_SPECFEM3D_SESAME)
+#
+#the modele create a python class for the mesh: ex. profile=mesh()
+#and it export the files of the mesh needed by the partitioner of SESAME
+#
+#############################################################################
+#OUTPUT
+#The default output are 11 ASCII files:
+#__________________________________________________________________________________________
+#mesh_name='mesh_file' -> the file that contains the connectity of the all mesh
+# format:
+# number of elements
+# id_elements id_node1 id_node2 id_node3 id_node4 id_node5 id_node6 id_node7 id_node8
+# .....
+#
+#__________________________________________________________________________________________
+##nodecoord_name='nodes_coords_file' -> the file that contains the coordinates of the nodes of the all mesh
+# format:
+# number of nodes
+# id_node x_coordinate y_coordinate z_coordinate
+# .....
+#
+#__________________________________________________________________________________________
+##material_name='materials_file' -> the file that contains the material flag of the elements
+# format:
+# id_element flag
+# .....
+#
+#__________________________________________________________________________________________
+##nummaterial_name='nummaterial_velocity_file' -> table of the material properties
+# format:
+# flag rho vp vs 0 0 #full definition of the properties, flag > 0
+# .....
+# flag 'tomography' file_name #for interpolation with tomography
+# .....
+# flag 'interface' file_name flag_for_the_gll_below_the_interface
+# flag_for_the_gll_above_the_interface #for interpolation with interface
+#__________________________________________________________________________________________
+##absname='absorbing_surface_file' -> this file contains all the face in all the absorbing boundaries
+##absname_local='absorbing_surface_file'+'_xmin' -> this file contains all the face in the
+# absorbing boundary defined by x=Xmin
+##absname_local='absorbing_surface_file'+'_xmax' -> this file contains all the face in the
+# absorbing boundary defined by x=Xmax
+##absname_local='absorbing_surface_file'+'_ymin' -> this file contains all the face in the
+# absorbing boundary defined by y=Ymin
+##absname_local='absorbing_surface_file'+'_ymax' -> this file contains all the face in the
+# absorbing boundary defined by y=Ymax
+##absname_local='absorbing_surface_file'+'_bottom' -> this file contains all the face in the
+# absorbing boundary defined by z=bottom
+# format:
+# number of faces
+# id_(element containg the face) id_node1_face id_node2_face id_node3_face id_node4_face
+# ....
+#
+#__________________________________________________________________________________________
+##freename='free_surface_file' -> file with the hex on the free surface (usually the topography)
+# format:
+# number of faces
+# id_(element containg the face) id_node1_face id_node2_face id_node3_face id_node4_face
+#
+#__________________________________________________________________________________________
+# it is possible save only one (or more) file singularly: for example if you want only the nodecoord_file
+# call the module mesh.nodescoord_write(full path name)
+#
+#############################################################################
+
+import cubit
+
+class mtools(object):
+ """docstring for ciao"""
+ def __init__(self,frequency,list_surf,list_vp):
+ super(mtools, self).__init__()
+ self.frequency = frequency
+ self.list_surf = list_surf
+ self.list_vp = list_vp
+ self.ngll=5
+ self.percent_gll=0.172
+ self.point_wavelength=5
+ def __repr__(self):
+ txt='Meshing for frequency up to '+str(self.frequency)+'Hz\n'
+ for surf,vp in zip(self.list_surf,self.list_vp):
+ txt=txt+'surface '+str(surf)+', vp ='+str(vp)+' -> size '+str(self.freq2meshsize(vp)[0])\
+ +' -> dt '+str(self.freq2meshsize(vp)[0])+'\n'
+ return txt
+ def freq2meshsize(self,vp):
+ velocity=vp*.5
+ self.size=(1/2.5)*velocity/self.frequency*(self.ngll-1)/self.point_wavelength
+ self.dt=.4*self.size/vp*self.percent_gll
+ return self.size,self.dt
+ def mesh_it(self):
+ for surf,vp in zip(self.list_surf,self.list_vp):
+ command = "surface "+str(surf)+" size "+str(self.freq2meshsize(vp)[0])
+ cubit.cmd(command)
+ command = "surface "+str(surf)+ 'scheme pave'
+ cubit.cmd(command)
+ command = "mesh surf "+str(surf)
+ cubit.cmd(command)
+
+class block_tools:
+ def __int__(self):
+ pass
+ def create_blocks(self,mesh_entity,list_entity=None,):
+ if mesh_entity =='surface':
+ txt=' face in surface '
+ elif mesh_entity == 'curve':
+ txt=' edge in curve '
+ elif mesh_entity == 'group':
+ txt=' face in group '
+ if list_entity:
+ if not isinstance(list_entity,list):
+ list_entity=[list_entity]
+ for entity in list_entity:
+ iblock=cubit.get_next_block_id()
+ command = "block "+str(iblock)+ txt +str(entity)
+ cubit.cmd(command)
+ def material_file(self,filename):
+ matfile=open(filename,'w')
+ material=[]
+ for record in matfile:
+ mat_name,vp_str=record.split()
+ vp=float(vp_str)
+ material.append([mat_name,vp])
+ self.material=dict(material)
+ def assign_block_material(self,id_block,mat_name,vp=None):
+ try:
+ material=self.material
+ except:
+ material=None
+ cubit.cmd('block '+str(id_block)+' attribute count 2')
+ cubit.cmd('block '+str(id_block)+' attribute index 1 '+str(id_block))
+ if material:
+ if material.has_key(mat_name):
+ cubit.cmd('block '+str(id_block)+' attribute index 2 '+str(material[mat_name]))
+ print 'block '+str(id_block)+' - material '+mat_name+' - vp '+str(material[mat_name])+' from database'
+ elif vp:
+ cubit.cmd('block '+str(id_block)+' attribute index 2 '+str(vp))
+ print 'block '+str(id_block)+' - material '+mat_name+' - vp '+str(vp)
+ else:
+ print 'assignment impossible: check if '+mat_name+' is in the database or specify vp'
+
+class mesh_tools(block_tools):
+ """Tools for the mesh
+ #########
+ dt,edge_dt,freq,edge_freq=seismic_resolution(edges,velocity,bins_d=None,bins_u=None,sidelist=None,ngll=5,np=8)
+ Given the velocity of a list of edges, seismic_resolution provides the minimum Dt
+ required for the stability condition (and the corrisponding edge).
+ Furthermore, given the number of gll point in the element (ngll) and the number
+ of GLL point for wavelength, it provide the maximum resolved frequency.
+ #########
+ length=edge_length(edge)
+ return the length of a edge
+ #########
+ edge_min,length=edge_min_length(surface)
+ given the cubit id of a surface, it return the edge with minimun length
+ #########
+ """
+ def __int__(self):
+ pass
+ def seismic_resolution(self,edges,velocity,bins_d=None,bins_u=None,sidelist=None):
+ """
+ dt,edge_dt,freq,edge_freq=seismic_resolution(edges,velocity,bins_d=None,bins_u=None,sidelist=None,ngll=5,np=8)
+ Given the velocity of a list of edges, seismic_resolution provides the minimum Dt
+ required for the stability condition (and the corrisponding edge).
+ Furthermore, given the number of gll point in the element (ngll) and the number
+ of GLL point for wavelength, it provide the maximum resolved frequency.
+ """
+ ratiostore=1e10
+ dtstore=1e10
+ edgedtstore=-1
+ edgeratiostore=-1
+ for edge in edges:
+ d=self.edge_length(edge)
+ ratio=(1/2.5)*velocity/d*(self.ngll-1)/self.point_wavelength
+ dt=.4*d/velocity*self.percent_gll
+ if dt<dtstore:
+ dtstore=dt
+ edgedtstore=edge
+ if ratio < ratiostore:
+ ratiostore=ratio
+ edgeratiostore=edge
+ try:
+ for bin_d,bin_u,side in zip(bins_d,bins_u,sidelist):
+ if ratio >= bin_d and ratio < bin_u:
+ command = "sideset "+str(side)+" edge "+str(edge)
+ cubit.cmd(command)
+ #print command
+ break
+ except:
+ pass
+ return dtstore,edgedtstore,ratiostore,edgeratiostore
+ def edge_length(self,edge):
+ """
+ length=edge_length(edge)
+ return the length of a edge
+ """
+ from math import sqrt
+ nodes=cubit.get_connectivity('Edge',edge)
+ x0,y0,z0=cubit.get_nodal_coordinates(nodes[0])
+ x1,y1,z1=cubit.get_nodal_coordinates(nodes[1])
+ d=sqrt((x1-x0)**2+(y1-y0)**2+(z1-z0)**2)
+ return d
+ def edge_min_length(self,surface):
+ """
+ edge_min,length=edge_min_length(surface)
+ given the cubit id of a surface, it return the edge with minimun length
+ """
+ from math import sqrt
+ self.dmin=99999
+ edge_store=0
+ command = "group 'list_edge' add edge in surf "+str(surface)
+ command = command.replace("["," ").replace("]"," ")
+ #print command
+ cubit.cmd(command)
+ group=cubit.get_id_from_name("list_edge")
+ edges=cubit.get_group_edges(group)
+ command = "delete group "+ str(group)
+ cubit.cmd(command)
+ for edge in edges:
+ d=self.edge_length(edge)
+ if d<dmin:
+ self.dmin=d
+ edge_store=edge
+ self.edgemin=edge_store
+ return self.edgemin,self.dmin
+ def normal_check(self,nodes,normal):
+ tres=.2
+ p0=cubit.get_nodal_coordinates(nodes[0])
+ p1=cubit.get_nodal_coordinates(nodes[1])
+ p2=cubit.get_nodal_coordinates(nodes[2])
+ a=[p1[0]-p0[0],p1[1]-p0[1],p1[2]-p0[2]]
+ b=[p2[0]-p1[0],p2[1]-p1[1],p2[2]-p1[2]]
+ axb=[a[1]*b[2] - a[2]*b[1], a[2]*b[0] - a[0]*b[2], a[0]*b[1] - a[1]*b[0]]
+ dot=0.0
+ for i in (0,1,2):
+ dot=dot+axb[i]*normal[i]
+ if dot > 0:
+ return nodes
+ elif dot < 0:
+ return nodes[0],nodes[3],nodes[2],nodes[1]
+ else:
+ print 'error: surface normal, dot=0', axb,normal,dot,p0,p1,p2
+ def mesh_analysis(self,frequency):
+ from sets import Set
+ cubit.cmd('set info off')
+ cubit.cmd('set echo off')
+ cubit.cmd('set journal off')
+ bins_d=[0.0001]+range(0,int(frequency)+1)+[1000]
+ bins_u=bins_d[1:]
+ dt=[]
+ ed_dt=[]
+ r=[]
+ ed_r=[]
+ nstart=cubit.get_next_sideset_id()
+ command = "del sideset all"
+ cubit.cmd(command)
+ for bin_d,bin_u in zip(bins_d,bins_u):
+ nsideset=cubit.get_next_sideset_id()
+ command='create sideset '+str(nsideset)
+ cubit.cmd(command)
+ command = "sideset "+str(nsideset)+ " name "+ "'ratio-["+str(bin_d)+"_"+str(bin_u)+"['"
+ cubit.cmd(command)
+ nend=cubit.get_next_sideset_id()
+ sidelist=range(nstart,nend)
+ for block in self.block_mat:
+ name=cubit.get_exodus_entity_name('block',block)
+ velocity=self.material[name][1]
+ if velocity > 0:
+ faces=cubit.get_block_faces(block)
+ edges=[]
+ for face in faces:
+ es=cubit.get_sub_elements("face", face, 1)
+ edges=edges+list(es)
+ edges=Set(edges)
+ dtstore,edgedtstore,ratiostore,edgeratiostore=self.seismic_resolution(edges,\
+ velocity,bins_d,bins_u,sidelist)
+ dt.append(dtstore)
+ ed_dt.append(edgedtstore)
+ r.append(ratiostore)
+ ed_r.append(edgeratiostore)
+ self.ddt=zip(ed_dt,dt)
+ self.dr=zip(ed_r,r)
+ def sorter(x, y):
+ return cmp(x[1],y[1])
+ self.ddt.sort(sorter)
+ self.dr.sort(sorter)
+ print self.ddt,self.dr
+ print 'Deltat minimum => edge:'+str(self.ddt[0][0])+' dt: '+str(self.ddt[0][1])
+ print 'minimum frequency resolved => edge:'+str(self.dr[0][0])+' frequency: '+str(self.dr[0][1])
+ return self.ddt[0],self.dr[0]
+
+class mesh(object,mesh_tools):
+ def __init__(self):
+ super(mesh, self).__init__()
+ self.mesh_name='mesh_file'
+ self.nodecoord_name='nodes_coords_file'
+ self.material_name='materials_file'
+ self.nummaterial_name='nummaterial_velocity_file'
+ self.absname='absorbing_surface_file'
+ self.freename='free_surface_file'
+ self.recname='STATIONS'
+ self.face='QUAD4'
+ self.face2='SHELL4'
+ self.hex='HEX8'
+ self.edge='BAR2'
+ self.topo='face_topo'
+ self.rec='receivers'
+ self.ngll=5
+ self.percent_gll=0.172
+ self.point_wavelength=5
+ self.block_definition()
+ cubit.cmd('compress')
+ def __repr__(self):
+ pass
+ def block_definition(self):
+ block_flag=[]
+ block_mat=[]
+ block_bc=[]
+ block_bc_flag=[]
+ material={}
+ bc={}
+ blocks=cubit.get_block_id_list()
+ for block in blocks:
+ name=cubit.get_exodus_entity_name('block',block)
+ type=cubit.get_block_element_type(block)
+ print block,name,blocks,type,self.hex,self.face
+ # block has hexahedral elements (HEX8)
+ if type == self.hex:
+ flag=None
+ vel=None
+ vs=None
+ rho=None
+ q=0
+ ani=0
+ # material domain id
+ if name.find("acoustic") >= 0 :
+ imaterial = 1
+ elif name.find("elastic") >= 0 :
+ imaterial = 2
+ elif name.find("poroelastic") >= 0 :
+ imaterial = 3
+ else :
+ imaterial = 0
+ print "block: ",name
+ print " could not find appropriate material for this block..."
+ print ""
+ break
+
+ nattrib=cubit.get_block_attribute_count(block)
+ if nattrib != 0:
+ # material flag:
+ # positive => material properties,
+ # negative => interface/tomography domain
+ flag=int(cubit.get_block_attribute_value(block,0))
+ if flag > 0 and nattrib >= 2:
+ # vp
+ vel=cubit.get_block_attribute_value(block,1)
+ if nattrib >= 3:
+ # vs
+ vs=cubit.get_block_attribute_value(block,2)
+ if nattrib >= 4:
+ #density
+ rho=cubit.get_block_attribute_value(block,3)
+ if nattrib >= 5:
+ #Q_flag
+ q=cubit.get_block_attribute_value(block,4)
+ # for q to be valid: it must be an integer flag between 1 and 13
+ # (see constants.h for IATTENUATION_SEDIMENT_40, etc. )
+ if q < 0 or q > 13:
+ print 'error, q flag invalid:', q
+ print ' check with constants.h for IATTENUATION flags'
+ break
+ if nattrib == 6:
+ #anisotropy_flag
+ ani=cubit.get_block_attribute_value(block,5)
+ elif flag < 0:
+ # velocity model
+ vel=name
+ attrib=cubit.get_block_attribute_value(block,1)
+ if attrib == 1:
+ kind='interface'
+ flag_down=cubit.get_block_attribute_value(block,2)
+ flag_up=cubit.get_block_attribute_value(block,3)
+ elif attrib == 2:
+ kind='tomography'
+ else:
+ flag=block
+ vel,vs,rho,q,ani=(name,0,0,0,0)
+ block_flag.append(int(flag))
+ block_mat.append(block)
+ if flag > 0:
+ par=tuple([imaterial,flag,vel,vs,rho,q,ani])
+ elif flag < 0:
+ if kind=='interface':
+ par=tuple([imaterial,flag,kind,name,flag_down,flag_up])
+ elif kind=='tomography':
+ par=tuple([imaterial,flag,kind,name])
+ elif flag==0:
+ par=tuple([imaterial,flag,name])
+ material[block]=par
+ elif (type == self.face) or (type == self.face2) :
+ # block has surface elements (QUAD4 or SHELL4)
+ block_bc_flag.append(4)
+ block_bc.append(block)
+ bc[block]=4 #face has connectivity = 4
+ if name == self.topo: topography_face=block
+ else:
+ # block elements differ from HEX8/QUAD4/SHELL4
+ print '****************************************'
+ print 'block not properly defined:'
+ print ' name:',name
+ print ' type:',type
+ print
+ print 'please check your block definitions!'
+ print
+ print 'only supported types are:'
+ print ' HEX8 for volumes'
+ print ' QUAD4 for surface'
+ print ' SHELL4 for surface'
+ print '****************************************'
+ continue
+
+ nsets=cubit.get_nodeset_id_list()
+ if len(nsets) == 0: self.receivers=None
+ for nset in nsets:
+ name=cubit.get_exodus_entity_name('nodeset',nset)
+ if name == self.rec:
+ self.receivers=nset
+ else:
+ print 'nodeset '+name+' not defined'
+ self.receivers=None
+ try:
+ self.block_mat=block_mat
+ self.block_flag=block_flag
+ self.block_bc=block_bc
+ self.block_bc_flag=block_bc_flag
+ self.material=material
+ self.bc=bc
+ self.topography=topography_face
+ except:
+ print '****************************************'
+ print 'sorry, no blocks or blocks not properly defined'
+ print block_mat
+ print block_flag
+ print block_bc
+ print block_bc_flag
+ print material
+ print bc
+ print topography
+ print '****************************************'
+ def mat_parameter(self,properties):
+ #note: material property acoustic/elastic/poroelastic are defined by the block's name
+ print "#material properties:"
+ print properties
+ imaterial=properties[0]
+ flag=properties[1]
+ if flag > 0:
+ vel=properties[2]
+ if properties[3] is None and type(vel) != str:
+ # velocity model scales with given vp value
+ if vel >= 30:
+ m2km=1000.
+ else:
+ m2km=1.
+ vp=vel/m2km
+ rho=(1.6612*vp-0.472*vp**2+0.0671*vp**3-0.0043*vp**4+0.000106*vp**4)*m2km
+ txt='%1i %3i %20f %20f %20f %1i %1i\n' % (properties[0],properties[1],rho,vel,vel/(3**.5),0,0)
+ elif type(vel) != str:
+ # velocity model given as vp,vs,rho,..
+ #format nummaterials file: #material_domain_id #material_id #rho #vp #vs #Q_flag #anisotropy_flag
+ txt='%1i %3i %20f %20f %20f %2i %2i\n' % (properties[0],properties[1],properties[4], \
+ properties[2],properties[3],properties[5],properties[6])
+ else:
+ txt='%1i %3i %s \n' % (properties[0],properties[1],properties[2])
+ elif flag < 0:
+ if properties[2] == 'tomography':
+ txt='%1i %3i %s %s\n' % (properties[0],properties[1],properties[2],properties[3])
+ elif properties[2] == 'interface':
+ txt='%1i %3i %s %s %1i %1i\n' % (properties[0],properties[1],properties[2],properties[3],\
+ properties[4],properties[5])
+ return txt
+ def nummaterial_write(self,nummaterial_name):
+ print 'Writing '+nummaterial_name+'.....'
+ nummaterial=open(nummaterial_name,'w')
+ for block in self.block_mat:
+ #name=cubit.get_exodus_entity_name('block',block)
+ nummaterial.write(self.mat_parameter(self.material[block]))
+ nummaterial.close()
+ print 'Ok'
+ def mesh_write(self,mesh_name):
+ meshfile=open(mesh_name,'w')
+ print 'Writing '+mesh_name+'.....'
+ num_elems=cubit.get_hex_count()
+ print num_elems
+ meshfile.write(str(num_elems)+'\n')
+ num_write=0
+ for block,flag in zip(self.block_mat,self.block_flag):
+ #print block,flag
+ hexes=cubit.get_block_hexes(block)
+ #print len(hexes)
+ for hexa in hexes:
+ #print hexa
+ nodes=cubit.get_connectivity('Hex',hexa)
+ #nodes=self.jac_check(nodes) #is it valid for 3D? TODO
+ txt=('%10i ')% hexa
+ txt=txt+('%10i %10i %10i %10i %10i %10i %10i %10i\n')% nodes[:]
+ meshfile.write(txt)
+ meshfile.close()
+ print 'Ok',str(num_elems)
+ def material_write(self,mat_name):
+ mat=open(mat_name,'w')
+ print 'Writing '+mat_name+'.....'
+ for block,flag in zip(self.block_mat,self.block_flag):
+ hexes=cubit.get_block_hexes(block)
+ for hexa in hexes:
+ mat.write(('%10i %10i\n') % (hexa,flag))
+ mat.close()
+ print 'Ok'
+ def nodescoord_write(self,nodecoord_name):
+ nodecoord=open(nodecoord_name,'w')
+ print 'Writing '+nodecoord_name+'.....'
+ node_list=cubit.parse_cubit_list('node','all')
+ num_nodes=len(node_list)
+ nodecoord.write('%10i\n' % num_nodes)
+ #
+ for node in node_list:
+ x,y,z=cubit.get_nodal_coordinates(node)
+ txt=('%10i %20f %20f %20f\n') % (node,x,y,z)
+ nodecoord.write(txt)
+ nodecoord.close()
+ print 'Ok'
+ def free_write(self,freename=None):
+ # free surface
+ cubit.cmd('set info off')
+ cubit.cmd('set echo off')
+ cubit.cmd('set journal off')
+ from sets import Set
+ normal=(0,0,1)
+ if not freename: freename=self.freename
+ # writes free surface file
+ print 'Writing '+freename+'.....'
+ freehex=open(freename,'w')
+ # searches block definition with name face_topo
+ for block,flag in zip(self.block_bc,self.block_bc_flag):
+ if block == self.topography:
+ name=cubit.get_exodus_entity_name('block',block)
+ print ' block name:',name,'id:',block
+ quads_all=cubit.get_block_faces(block)
+ print ' face = ',len(quads_all)
+ dic_quads_all=dict(zip(quads_all,quads_all))
+ freehex.write('%10i\n' % len(quads_all))
+ list_hex=cubit.parse_cubit_list('hex','all')
+ for h in list_hex:
+ faces=cubit.get_sub_elements('hex',h,2)
+ for f in faces:
+ if dic_quads_all.has_key(f):
+ #print f
+ nodes=cubit.get_connectivity('Face',f)
+ nodes_ok=self.normal_check(nodes,normal)
+ txt='%10i %10i %10i %10i %10i\n' % (h,nodes_ok[0],\
+ nodes_ok[1],nodes_ok[2],nodes_ok[3])
+ freehex.write(txt)
+ freehex.close()
+ print 'Ok'
+ cubit.cmd('set info on')
+ cubit.cmd('set echo on')
+ def abs_write(self,absname=None):
+ # absorbing boundaries
+ import re
+ cubit.cmd('set info off')
+ cubit.cmd('set echo off')
+ cubit.cmd('set journal off')
+ from sets import Set
+ if not absname: absname=self.absname
+ #
+ # loops through all block definitions
+ list_hex=cubit.parse_cubit_list('hex','all')
+ for block,flag in zip(self.block_bc,self.block_bc_flag):
+ if block != self.topography:
+ name=cubit.get_exodus_entity_name('block',block)
+ print name,block
+ absflag=False
+ if re.search('xmin',name):
+ filename=absname+'_xmin'
+ normal=(-1,0,0)
+ elif re.search('xmax',name):
+ filename=absname+'_xmax'
+ normal=(1,0,0)
+ elif re.search('ymin',name):
+ filename=absname+'_ymin'
+ normal=(0,-1,0)
+ elif re.search('ymax',name):
+ filename=absname+'_ymax'
+ normal=(0,1,0)
+ elif re.search('bottom',name):
+ filename=absname+'_bottom'
+ normal=(0,0,-1)
+ elif re.search('abs',name):
+ print " ...face_abs - not used so far..."
+ continue
+ else:
+ continue
+ # opens file
+ print 'Writing '+filename+'.....'
+ abshex_local=open(filename,'w')
+ # gets face elements
+ quads_all=cubit.get_block_faces(block)
+ dic_quads_all=dict(zip(quads_all,quads_all))
+ abshex_local.write('%10i\n' % len(quads_all))
+ #command = "group 'list_hex' add hex in face "+str(quads_all)
+ #command = command.replace("["," ").replace("]"," ").replace("("," ").replace(")"," ")
+ #cubit.cmd(command)
+ #group=cubit.get_id_from_name("list_hex")
+ #list_hex=cubit.get_group_hexes(group)
+ #command = "delete group "+ str(group)
+ #cubit.cmd(command)
+ for h in list_hex:
+ faces=cubit.get_sub_elements('hex',h,2)
+ for f in faces:
+ if dic_quads_all.has_key(f):
+ nodes=cubit.get_connectivity('Face',f)
+ if not absflag:
+ # checks with specified normal
+ nodes_ok=self.normal_check(nodes,normal)
+ txt='%10i %10i %10i %10i %10i\n' % (h,nodes_ok[0],\
+ nodes_ok[1],nodes_ok[2],nodes_ok[3])
+ else:
+ txt='%10i %10i %10i %10i %10i\n' % (h,nodes[0],\
+ nodes[1],nodes[2],nodes[3])
+ abshex_local.write(txt)
+ # closes file
+ abshex_local.close()
+ print 'Ok'
+ cubit.cmd('set info on')
+ cubit.cmd('set echo on')
+ def surface_write(self,pathdir=None):
+ # optional surfaces, e.g. moho_surface
+ # should be created like e.g.:
+ # > block 10 face in surface 2
+ # > block 10 name 'moho_surface'
+ import re
+ from sets import Set
+ for block in self.block_bc :
+ if block != self.topography:
+ name=cubit.get_exodus_entity_name('block',block)
+ # skips block names like face_abs**, face_topo**
+ if re.search('abs',name):
+ continue
+ elif re.search('topo',name):
+ continue
+ elif re.search('surface',name):
+ filename=pathdir+name+'_file'
+ else:
+ continue
+ # gets face elements
+ print ' surface block name: ',name,'id: ',block
+ quads_all=cubit.get_block_faces(block)
+ print ' face = ',len(quads_all)
+ if len(quads_all) == 0 :
+ continue
+ # writes out surface infos to file
+ print 'Writing '+filename+'.....'
+ surfhex_local=open(filename,'w')
+ dic_quads_all=dict(zip(quads_all,quads_all))
+ # writes number of surface elements
+ surfhex_local.write('%10i\n' % len(quads_all))
+ # writes out element node ids
+ list_hex=cubit.parse_cubit_list('hex','all')
+ for h in list_hex:
+ faces=cubit.get_sub_elements('hex',h,2)
+ for f in faces:
+ if dic_quads_all.has_key(f):
+ nodes=cubit.get_connectivity('Face',f)
+ txt='%10i %10i %10i %10i %10i\n' % (h,nodes[0],\
+ nodes[1],nodes[2],nodes[3])
+ surfhex_local.write(txt)
+ # closes file
+ surfhex_local.close()
+ print 'Ok'
+ def rec_write(self,recname):
+ print 'Writing '+self.recname+'.....'
+ recfile=open(self.recname,'w')
+ nodes=cubit.get_nodeset_nodes(self.receivers)
+ for i,n in enumerate(nodes):
+ x,y,z=cubit.get_nodal_coordinates(n)
+ recfile.write('ST%i XX %20f %20f 0.0 0.0 \n' % (i,x,z))
+ recfile.close()
+ print 'Ok'
+ def write(self,path=''):
+ cubit.cmd('set info off')
+ cubit.cmd('set echo off')
+ cubit.cmd('set journal off')
+ if len(path) != 0:
+ if path[-1] != '/': path=path+'/'
+ # mesh file
+ self.mesh_write(path+self.mesh_name)
+ # mesh material
+ self.material_write(path+self.material_name)
+ # mesh coordinates
+ self.nodescoord_write(path+self.nodecoord_name)
+ # material definitions
+ self.nummaterial_write(path+self.nummaterial_name)
+ # free surface: face_top
+ self.free_write(path+self.freename)
+ # absorbing surfaces: abs_***
+ self.abs_write(path+self.absname)
+ # any other surfaces: ***surface***
+ self.surface_write(path)
+ # receivers
+ if self.receivers: self.rec_write(path+self.recname)
+ cubit.cmd('set info on')
+ cubit.cmd('set echo on')
+
+def export2SESAME(path_exporting_mesh_SPECFEM3D_SESAME):
+ cubit.cmd('set info on')
+ cubit.cmd('set echo on')
+ sem_mesh=mesh()
+ sem_mesh.write(path=path_exporting_mesh_SPECFEM3D_SESAME)
+
+
+if __name__ == '__main__':
+ path='MESH/'
+ export2SESAME(path)
+
+# call by:
+# import cubit2specfem3d
+# cubit2specfem3d.export2SESAME('MESH')
Added: seismo/3D/FAULT_SOURCE/branches/new_fault_db/CUBIT/fault_elements_1.dat
===================================================================
--- seismo/3D/FAULT_SOURCE/branches/new_fault_db/CUBIT/fault_elements_1.dat (rev 0)
+++ seismo/3D/FAULT_SOURCE/branches/new_fault_db/CUBIT/fault_elements_1.dat 2011-05-09 15:41:36 UTC (rev 18335)
@@ -0,0 +1,127 @@
+ 63 63
+ 36 52 44 147 170
+ 37 51 52 170 169
+ 38 50 51 169 168
+ 39 49 50 168 167
+ 40 48 49 167 166
+ 41 47 48 166 165
+ 42 46 47 165 164
+ 43 45 46 164 163
+ 44 43 45 163 162
+ 146 170 147 146 178
+ 147 169 170 178 177
+ 148 168 169 177 176
+ 149 167 168 176 175
+ 150 166 167 175 174
+ 151 165 166 174 173
+ 152 164 165 173 172
+ 153 163 164 172 171
+ 154 162 163 171 161
+ 256 178 146 145 186
+ 257 177 178 186 185
+ 258 176 177 185 184
+ 259 175 176 184 183
+ 260 174 175 183 182
+ 261 173 174 182 181
+ 262 172 173 181 180
+ 263 171 172 180 179
+ 264 161 171 179 160
+ 366 186 145 144 194
+ 367 185 186 194 193
+ 368 184 185 193 192
+ 369 183 184 192 191
+ 370 182 183 191 190
+ 371 181 182 190 189
+ 372 180 181 189 188
+ 373 179 180 188 187
+ 374 160 179 187 159
+ 476 194 144 143 202
+ 477 193 194 202 201
+ 478 192 193 201 200
+ 479 191 192 200 199
+ 480 190 191 199 198
+ 481 189 190 198 197
+ 482 188 189 197 196
+ 483 187 188 196 195
+ 484 159 187 195 158
+ 586 202 143 142 210
+ 587 201 202 210 209
+ 588 200 201 209 208
+ 589 199 200 208 207
+ 590 198 199 207 206
+ 591 197 198 206 205
+ 592 196 197 205 204
+ 593 195 196 204 203
+ 594 158 195 203 157
+ 696 210 142 141 156
+ 697 209 210 156 155
+ 698 208 209 155 154
+ 699 207 208 154 153
+ 700 206 207 153 152
+ 701 205 206 152 151
+ 702 204 205 151 150
+ 703 203 204 150 149
+ 704 157 203 149 148
+ 25 60 43 162 226
+ 26 59 60 226 225
+ 27 58 59 225 224
+ 28 57 58 224 223
+ 29 56 57 223 222
+ 30 55 56 222 221
+ 31 54 55 221 220
+ 32 53 54 220 219
+ 33 44 53 219 147
+ 135 226 162 161 234
+ 136 225 226 234 233
+ 137 224 225 233 232
+ 138 223 224 232 231
+ 139 222 223 231 230
+ 140 221 222 230 229
+ 141 220 221 229 228
+ 142 219 220 228 227
+ 143 147 219 227 146
+ 245 234 161 160 242
+ 246 233 234 242 241
+ 247 232 233 241 240
+ 248 231 232 240 239
+ 249 230 231 239 238
+ 250 229 230 238 237
+ 251 228 229 237 236
+ 252 227 228 236 235
+ 253 146 227 235 145
+ 355 242 160 159 250
+ 356 241 242 250 249
+ 357 240 241 249 248
+ 358 239 240 248 247
+ 359 238 239 247 246
+ 360 237 238 246 245
+ 361 236 237 245 244
+ 362 235 236 244 243
+ 363 145 235 243 144
+ 465 250 159 158 258
+ 466 249 250 258 257
+ 467 248 249 257 256
+ 468 247 248 256 255
+ 469 246 247 255 254
+ 470 245 246 254 253
+ 471 244 245 253 252
+ 472 243 244 252 251
+ 473 144 243 251 143
+ 575 258 158 157 266
+ 576 257 258 266 265
+ 577 256 257 265 264
+ 578 255 256 264 263
+ 579 254 255 263 262
+ 580 253 254 262 261
+ 581 252 253 261 260
+ 582 251 252 260 259
+ 583 143 251 259 142
+ 685 266 157 148 211
+ 686 265 266 211 212
+ 687 264 265 212 213
+ 688 263 264 213 214
+ 689 262 263 214 215
+ 690 261 262 215 216
+ 691 260 261 216 217
+ 692 259 260 217 218
+ 693 142 259 218 141
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-08 08:12:31 UTC (rev 18334)
+++ seismo/3D/FAULT_SOURCE/branches/new_fault_db/decompose_mesh_SCOTCH/fault_scotch.f90 2011-05-09 15:41:36 UTC (rev 18335)
@@ -53,21 +53,25 @@
integer,intent(in) :: ifault
character(len=5) :: NTchar
- integer :: e,ier
-
+ integer :: e,ier,nspec_side1,nspec_side2
+
write(NTchar,'(I5)') ifault
NTchar = adjustl(NTchar)
- open(101,file='../DATA/FAULT/fault_elements_'//NTCHAR//'.dat',&
+ open(101,file='../DATA/FAULT/fault_file_'//NTCHAR//'.dat',&
status='old',action='read',iostat=ier)
if( ier == 0 ) then
+ read(101,*) nspec_side1,nspec_side2
+ if (nspec_side1 =/ nspec_side2) stop 'Number of fault nodes at do not match.'
+ f%nspec = nspec_side1
read(101,*) f%nspec
allocate(f%ispec1(f%nspec))
allocate(f%ispec2(f%nspec))
allocate(f%inodes1(4,f%nspec))
allocate(f%inodes2(4,f%nspec))
! format:
+ ! number of elements side at side 1 and side 2
! #id_(element containing the face) #id_node1_face .. #id_node4_face
! First for all faces on side 1, then side 2
! Note: element ids start at 1, not 0. Check that in cubit2specfem3d.py
@@ -84,7 +88,7 @@
! read(101,*) f%ispec1(e),f%ispec2(e),f%iface1(e),f%iface2(e)
!enddo
else
- write(6,*) 'Fatal error: file ../DATA/FAULT/fault_elements_'//NTCHAR//'.dat not found'
+ write(6,*) 'Fatal error: file ../DATA/FAULT/fault_file_'//NTCHAR//'.dat not found'
write(6,*) 'Abort'
stop
endif
@@ -94,7 +98,6 @@
! ---------------------------------------------------------------------------------------------------
-!
subroutine close_faults(nodes_coords,elmnts,nelmnts,nnodes,esize)
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-08 08:12:31 UTC (rev 18334)
+++ seismo/3D/FAULT_SOURCE/branches/new_fault_db/src/fault_object.f90 2011-05-09 15:41:36 UTC (rev 18335)
@@ -268,7 +268,7 @@
type(fault_db_type), intent(in) :: fdb
integer, intent(in) :: nspec ! number of spectral elements in each block
- if (fdb(iflt)%eta > 0.0_CUSTOM_REAL) then
+ if (fdb%eta > 0.0_CUSTOM_REAL) then
if (.not.allocated(Kelvin_Voigt_eta)) then
allocate(Kelvin_Voigt_eta(nspec))
Kelvin_Voigt_eta(:) = 0.0_CUSTOM_REAL
More information about the CIG-COMMITS
mailing list