From d2a344be186677d2394c955d9c79354c8ff4bcca Mon Sep 17 00:00:00 2001 From: Daniel Peter Date: Tue, 29 Oct 2019 13:00:06 +0300 Subject: [PATCH] updates python print statements for future use with python3 --- CUBIT_GEOCUBIT/GEOCUBIT.py | 24 +- .../geocubitlib/LatLongUTMconversion.py | 7 +- .../geocubitlib/boundary_definition.py | 100 ++-- CUBIT_GEOCUBIT/geocubitlib/cubit2specfem3d.py | 255 ++++----- CUBIT_GEOCUBIT/geocubitlib/dev/cmpl.py | 11 +- CUBIT_GEOCUBIT/geocubitlib/exportlib.py | 191 +++---- CUBIT_GEOCUBIT/geocubitlib/hex_metric.py | 146 ++--- CUBIT_GEOCUBIT/geocubitlib/local_volume.py | 28 +- CUBIT_GEOCUBIT/geocubitlib/menu.py | 42 +- CUBIT_GEOCUBIT/geocubitlib/mesh_volume.py | 35 +- CUBIT_GEOCUBIT/geocubitlib/quality_log.py | 8 +- .../geocubitlib/read_parameter_cfg.py | 14 +- CUBIT_GEOCUBIT/geocubitlib/sea.py | 1 + CUBIT_GEOCUBIT/geocubitlib/start.py | 26 +- CUBIT_GEOCUBIT/geocubitlib/surfaces.py | 10 +- CUBIT_GEOCUBIT/geocubitlib/utilities.py | 28 +- CUBIT_GEOCUBIT/geocubitlib/volumes.py | 62 +-- .../block_mesh-anisotropic.py | 1 + .../block_mesh-anisotropic.py | 1 + .../block_mesh-anisotropic.py | 1 + .../block_mesh-anisotropic.py | 1 + .../block_mesh-anisotropic.py | 1 + .../block_mesh-anisotropic.py | 1 + EXAMPLES/Gmsh_simple_lddrk/Gmsh2specfem.py | 510 ------------------ EXAMPLES/Gmsh_simple_lddrk/README | 4 +- .../mesh_example_with_gmsh.py | 50 +- .../block_mesh.py | 3 +- EXAMPLES/Mount_StHelens/mesh_mount.py | 92 ++-- EXAMPLES/Mount_StHelens/mesh_mount_stl.py | 55 +- EXAMPLES/Mount_StHelens/read_topo.py | 33 +- .../splay_faults/splay_faults.py | 1 + EXAMPLES/fault_examples/tpv102/TPV102.py | 1 + EXAMPLES/fault_examples/tpv103/TPV103.py | 1 + EXAMPLES/fault_examples/tpv15/TPV15.py | 1 + EXAMPLES/fault_examples/tpv16/TPV16.py | 1 + EXAMPLES/fault_examples/tpv5/TPV5.py | 2 + .../block_mesh-anisotropic.py | 1 + .../create_adjoint_sources.sh | 3 + .../create_alpha_kernel_vtk.sh | 3 + .../create_mesh.py | 62 +-- .../make_mesh.sh | 8 +- .../block_mesh-anisotropic.py | 9 +- .../2lay_mesh_boundary_fig8-nodoubling.py | 4 +- .../2lay_mesh_boundary_fig8.py | 4 +- .../DATA/py/forcesolution_generation.py | 2 + .../cavity/DATA/py/stations_generation.py | 2 + EXAMPLES/tomographic_model/tomoblock_mesh.py | 22 +- .../make_mesh_waterlayer-default.py | 4 +- .../make_mesh_waterlayer_verticalboundary.py | 1 + .../waterlayer_mesh_boundary_fig8.py | 2 + .../waterlayer_mesh_boundary_vertical.py | 1 + .../waterlayer_only-nodoubling.py | 3 +- .../input_output/mesh_tools_mod.f90 | 4 +- ...bGmsh2Specfem_convert_Gmsh_to_Specfem3D.py | 74 +-- .../GMT/grenoble_basin_GMT_movies/mkmpeg4.py | 21 +- .../Paraview/paraviewpython-example.py | 21 +- utils/Visualization/plot_beachball.py | 1 + ...script_convert_signals_to_single_binary.py | 2 + .../test_script_open_binarized_file.py | 2 + .../Present_specfem3D_results.py | 28 +- .../curved_fault_mesh/create_box_mesh.py | 48 +- .../create_semisphere_mesh.py | 34 +- .../python_scripts/Dimension.py | 4 +- .../python_scripts/Smooth_slab.py | 4 +- .../python_scripts/playback_jou.py | 8 +- .../python_scripts/slab_interface_netsurf.py | 8 +- .../dipping_fault_planar_dip10_kink.py | 32 +- .../stepover/change_to_excel_format.py | 5 +- .../stepover/check_jump_or_not.py | 22 +- .../stepover/edit_parfile_fault.py | 10 +- .../stepover/set_uniform_stress.py | 4 +- .../subductionzone_mesh/exportmesh.py | 2 + .../process_slab_rotate.py | 2 + 73 files changed, 930 insertions(+), 1285 deletions(-) delete mode 100755 EXAMPLES/Gmsh_simple_lddrk/Gmsh2specfem.py diff --git a/CUBIT_GEOCUBIT/GEOCUBIT.py b/CUBIT_GEOCUBIT/GEOCUBIT.py index d21bc30e1..bf471f167 100755 --- a/CUBIT_GEOCUBIT/GEOCUBIT.py +++ b/CUBIT_GEOCUBIT/GEOCUBIT.py @@ -30,6 +30,24 @@ - numpy 1.0+ - http://downloads.sourceforge.net/numpy """ +from __future__ import print_function + +import sys + +print("GEOCUBIT") + +# version info +python_major_version = sys.version_info[0] +python_minor_version = sys.version_info[1] +print("Python version: ","{}.{}".format(python_major_version,python_minor_version)) + +# in case importing menu fails due to import utilities errors to find, +# this will add the geocubitlib/ folder to the sys.path: +import geocubitlib +sys.path.append(geocubitlib.__path__[0]) +#print(sys.path) +print('') + import geocubitlib.menu as menu import geocubitlib.start as start @@ -41,7 +59,7 @@ if __name__ == '__main__': - print 'VERSION 4.1' + print('VERSION 4.1') # GEOMETRY if menu.build_surface: @@ -81,8 +99,8 @@ if menu.export2SPECFEM3D and not menu.collect: from geocubitlib.exportlib import e2SEM - print menu.cubfiles - print 'hex27 ', menu.hex27 + print(menu.cubfiles) + print('hex27 ', menu.hex27) e2SEM(files=menu.cubfiles, listblock=menu.listblock, listflag=menu.listflag, diff --git a/CUBIT_GEOCUBIT/geocubitlib/LatLongUTMconversion.py b/CUBIT_GEOCUBIT/geocubitlib/LatLongUTMconversion.py index 19be3c38e..48e32e8ab 100755 --- a/CUBIT_GEOCUBIT/geocubitlib/LatLongUTMconversion.py +++ b/CUBIT_GEOCUBIT/geocubitlib/LatLongUTMconversion.py @@ -1,6 +1,7 @@ #!/usr/bin/env python - +# # Lat Long - UTM, UTM - Lat Long conversions +from __future__ import print_function from math import pi, sin, cos, tan, sqrt @@ -246,6 +247,6 @@ def UTMtoLL(ReferenceEllipsoid, northing, easting, zone): if __name__ == '__main__': (z, e, n) = LLtoUTM(23, 45.00, -75.00) - print z, e, n + print(z, e, n) (lat, lon) = UTMtoLL(23, n, e, z) - print lat, lon + print(lat, lon) diff --git a/CUBIT_GEOCUBIT/geocubitlib/boundary_definition.py b/CUBIT_GEOCUBIT/geocubitlib/boundary_definition.py index bec45b5fc..805f28976 100755 --- a/CUBIT_GEOCUBIT/geocubitlib/boundary_definition.py +++ b/CUBIT_GEOCUBIT/geocubitlib/boundary_definition.py @@ -22,6 +22,8 @@ # 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. # # # ############################################################################# +from __future__ import print_function + try: import start as start cubit = start.start_cubit() @@ -29,7 +31,7 @@ try: import cubit except: - print 'error importing cubit, check if cubit is installed' + print('error importing cubit, check if cubit is installed') pass from utilities import list2str @@ -112,15 +114,15 @@ def lateral_boundary_are_absorbing(iproc=0, cpuxmin=0, cpuxmax=1, # if iproc in iproc_xmin: abs_xmin = xmin - print 'proc ', iproc, ' has absorbing boundary xmin' + print('proc ', iproc, ' has absorbing boundary xmin') if iproc in iproc_ymin: - print 'proc ', iproc, ' has absorbing boundary ymin' + print('proc ', iproc, ' has absorbing boundary ymin') abs_ymin = ymin if iproc in iproc_xmax: - print 'proc ', iproc, ' has absorbing boundary xmax' + print('proc ', iproc, ' has absorbing boundary xmax') abs_xmax = xmax if iproc in iproc_ymax: - print 'proc ', iproc, ' has absorbing boundary ymax' + print('proc ', iproc, ' has absorbing boundary ymax') abs_ymax = ymax return abs_xmin, abs_xmax, abs_ymin, abs_ymax @@ -279,9 +281,9 @@ def build_block(vol_list, name, id_0=1, top_surf=None, optionsea=False): seathres = False # - print 'build blocks' + print('build blocks') block_list = cubit.get_block_id_list() - print block_list, vol_list + print(block_list, vol_list) if len(block_list) > 0: id_block = max(max(block_list), 2) + id_0 else: @@ -318,11 +320,11 @@ def build_block(vol_list, name, id_0=1, top_surf=None, optionsea=False): if version_cubit >= 15: command = 'block ' + str(id_block) + ' hex in vol ' + \ str(v) - print command + print(command) else: command = 'block ' + str(id_block) + ' hex in vol ' + \ str(v) + ' except hex in vol ' + str(list(v_other)) - print command + print(command) command = command.replace("[", " ").replace("]", " ") cubit.cmd(command) command = "block " + str(id_block) + " name '" + n + "'" @@ -390,76 +392,76 @@ def define_bc(*args, **keys): cpux=cpux, cpuy=cpuy) id_0 = cubit.get_next_block_id() v_list, name_list = define_block() - print 'define block', v_list, name_list + print('define block', v_list, name_list) build_block(v_list, name_list, id_0, top_surf, optionsea=optionsea) # entities entities = ['face'] if type(args) == list: if len(args) > 0: entities = args[0] - print entities + print(entities) for entity in entities: - print "##entity: " + str(entity) + print("##entity: " + str(entity)) # block for free surface (w/ topography) - # print '## topo surface block: ' + str(topo) + # print('## topo surface block: ' + str(topo)) if len(top_surf) == 0: - print "" - print "no topo surface found,\ - please create block face_topo manually..." - print "" + print("") + print("no topo surface found,\ + please create block face_topo manually...") + print("") else: build_block_side(top_surf, entity + '_topo', obj=entity, id_0=1001) # model has parallel sides (e.g. a block model ) # xmin - blocks if len(xmin) == 0: - print "" - print "0 abs_xmin surface found, please create block manually" - print "" + print("") + print("0 abs_xmin surface found, please create block manually") + print("") else: build_block_side(xmin, entity + '_abs_xmin', obj=entity, id_0=1003) # xmax - blocks if len(xmax) == 0: - print "" - print "0 abs_xmax surface found, please create block manually" - print "" + print("") + print("0 abs_xmax surface found, please create block manually") + print("") else: build_block_side(xmax, entity + '_abs_xmax', obj=entity, id_0=1005) # ymin - blocks if len(ymin) == 0: - print "" - print "0 abs_xmin surface found, please create block manually" - print "" + print("") + print("0 abs_xmin surface found, please create block manually") + print("") else: build_block_side(ymin, entity + '_abs_ymin', obj=entity, id_0=1004) # ymax - blocks if len(ymax) == 0: - print "" - print "0 abs_ymax surface found, please create block manually" - print "" + print("") + print("0 abs_ymax surface found, please create block manually") + print("") else: build_block_side(ymax, entity + '_abs_ymax', obj=entity, id_0=1006) # bottom - blocks if len(bottom_surf) == 0: - print "" - print "0 abs_bottom surf found, please create block manually" - print "" + print("") + print("0 abs_bottom surf found, please create block manually") + print("") else: build_block_side(bottom_surf, entity + '_abs_bottom', obj=entity, id_0=1002) elif closed: - print "##closed region not ready" + print("##closed region not ready") # surf = define_absorbing_surf_sphere() # v_list, name_list = define_block() # build_block(v_list, name_list, id_0) # # entities # entities = args[0] # id_side = 1001 - # print entities + # print(entities) # for entity in entities: # build_block_side(surf, entity + '_closedvol', # obj=entity, id_0=id_side) @@ -528,7 +530,7 @@ def get_ordered_node_surf(lsurface, icurve): if k != 0: cubit.cmd('del group sl') else: - print 'initializing group sl' + print('initializing group sl') cubit.cmd("group 'sl' add node in surf " + lsurf) group1 = cubit.get_id_from_name("sl") nodes_ls = list(cubit.get_group_nodes(group1)) @@ -540,7 +542,7 @@ def get_ordered_node_surf(lsurface, icurve): if k != 0: cubit.cmd('del group n1') else: - print 'initializing group n1' + print('initializing group n1') cubit.cmd("group 'n1' add node in curve " + icurvestr) x = cubit.get_bounding_box('curve', icurve) if x[2] > x[5]: @@ -600,7 +602,7 @@ def get_ordered_node_surf(lsurface, icurve): if k != 0: cubit.cmd('del group curve_vertical') else: - print 'initializing group curve_vertical' + print('initializing group curve_vertical') cubit.cmd("group 'curve_vertical' add node in curve " + kcurve) group1 = cubit.get_id_from_name('curve_vertical') nodes_curve = list(cubit.get_group_nodes(group1)) @@ -632,14 +634,14 @@ def check_bc(iproc, xmin, xmax, ymin, ymax, curve_bottom_xmin, curve_bottom_ymin, curve_bottom_xmax, \ curve_bottom_ymax = \ extract_bottom_curves(surf_xmin, surf_ymin, surf_xmax, surf_ymax) - # print 'absorbing surfaces: ', absorbing_surf - print 'absorbing surfaces xmin : ', abs_xmin - print 'absorbing surfaces xmax : ', abs_xmax - print 'absorbing surfaces ymin : ', abs_ymin - print 'absorbing surfaces ymax : ', abs_ymax - print 'absorbing surfaces top : ', top_surf - print 'absorbing surfaces bottom : ', bottom_surf - # print 'bottom curves: ', surf_xmin, surf_ymin, surf_xmax, surf_ymax + # print('absorbing surfaces: ', absorbing_surf) + print('absorbing surfaces xmin : ', abs_xmin) + print('absorbing surfaces xmax : ', abs_xmax) + print('absorbing surfaces ymin : ', abs_ymin) + print('absorbing surfaces ymax : ', abs_ymax) + print('absorbing surfaces top : ', top_surf) + print('absorbing surfaces bottom : ', bottom_surf) + # print('bottom curves: ', surf_xmin, surf_ymin, surf_xmax, surf_ymax) # # # @@ -717,10 +719,10 @@ def check_bc(iproc, xmin, xmax, ymin, ymax, boundary['node_curve_xmaxymin'] = c_xmaxymin boundary['node_curve_xmaxymax'] = c_xmaxymax - # print boundary['node_curve_xminymin'] - # print boundary['node_curve_xminymax'] - # print boundary['node_curve_xmaxymin'] - # print boundary['node_curve_xmaxymax'] + # print(boundary['node_curve_xminymin']) + # print(boundary['node_curve_xminymax']) + # print(boundary['node_curve_xmaxymin']) + # print(boundary['node_curve_xmaxymax']) entities = ['face'] # for entity in entities: diff --git a/CUBIT_GEOCUBIT/geocubitlib/cubit2specfem3d.py b/CUBIT_GEOCUBIT/geocubitlib/cubit2specfem3d.py index 858cb73e7..6df84df88 100755 --- a/CUBIT_GEOCUBIT/geocubitlib/cubit2specfem3d.py +++ b/CUBIT_GEOCUBIT/geocubitlib/cubit2specfem3d.py @@ -127,6 +127,7 @@ # module mesh.nodescoord_write(full path name) # ############################################################################# +from __future__ import print_function try: import start as start @@ -135,7 +136,7 @@ try: import cubit except: - print 'error importing cubit, check if cubit is installed' + print('error importing cubit, check if cubit is installed') pass from utilities import get_cubit_version @@ -219,16 +220,16 @@ def assign_block_material(self, id_block, mat_name, vp=None): if mat_name in material.keys(): 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' + 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) + 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' + print('assignment impossible: check if ' + mat_name +\ + ' is in the database or specify vp') class mesh_tools(block_tools): @@ -345,7 +346,7 @@ def normal_check(self, nodes, normal): 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 + print('error: surface normal, dot=0', axb, normal, dot, p0, p1, p2) def mesh_analysis(self, frequency): from sets import Set @@ -396,11 +397,11 @@ def sorter(x, y): 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]) + 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] @@ -437,8 +438,8 @@ def __init__(self, hex27=False, cpml=False, cpml_size=False, self.hex = 'HEX' if version_cubit <= 13: if hex27: - print "ATTENTION **********************\n\nCubit <=\ - 12.2 doesn't support HEX27\nassuming HEX8 .....\n\n" + print("ATTENTION **********************\n\nCubit <=\ + 12.2 doesn't support HEX27\nassuming HEX8 .....\n\n") self.hex27 = False else: self.hex27 = hex27 @@ -453,7 +454,7 @@ def __init__(self, hex27=False, cpml=False, cpml_size=False, if cpml_size: self.size = cpml_size else: - print 'please specify cmpl size if you want to use cpml' + print('please specify cmpl size if you want to use cpml') self.top_absorbing = top_absorbing if hex27: cubit.cmd('block all except block 1001 1002 1003 1004 1005 1006 \ @@ -490,7 +491,7 @@ def block_definition(self): for block in blocks: name = cubit.get_exodus_entity_name('block', block) ty = cubit.get_block_element_type(block) - # print block,blocks,ty,self.hex,self.face + # print(block,blocks,ty,self.hex,self.face) if self.hex in ty: flag = None vp = None @@ -549,7 +550,7 @@ def block_definition(self): block, 5) # for q to be valid: it must be positive if qk < 0 or qmu < 0: - print 'error, Q value invalid:',qk, qmu + print('error, Q value invalid:',qk, qmu) break if nattrib == 7: # anisotropy_flag @@ -570,7 +571,7 @@ def block_definition(self): kind = 'tomography' elif nattrib == 1: flag = cubit.get_block_attribute_value(block, 0) - # print 'only 1 attribute ', name,block,flag + # print('only 1 attribute ', name,block,flag) vp, vs, rho, qk, qmu, ani = (0, 0, 0, 9999., 9999., 0) else: flag = block @@ -600,18 +601,18 @@ def block_definition(self): pass else: # block elements differ from HEX8/QUAD4/SHELL4 - print '****************************************' - print 'block not properly defined:' - print ' name:', name - print ' type:', ty - print - print 'please check your block definitions!' - print - print 'only supported types are:' - print ' HEX/HEX8/HEX27 for volumes' - print ' QUAD4 for surface' - print ' SHELL4/SHELL9 for surface' - print '****************************************' + print('****************************************') + print('block not properly defined:') + print(' name:', name) + print(' type:', ty) + print('') + print('please check your block definitions!') + print('') + print('only supported types are:') + print(' HEX/HEX8/HEX27 for volumes') + print(' QUAD4 for surface') + print(' SHELL4/SHELL9 for surface') + print('****************************************') continue return None, None, None, None, None, None, None, None nsets = cubit.get_nodeset_id_list() @@ -622,7 +623,7 @@ def block_definition(self): if name == self.rec: self.receivers = nset else: - print 'nodeset ' + name + ' not defined' + print('nodeset ' + name + ' not defined') self.receivers = None try: @@ -632,26 +633,26 @@ def block_definition(self): self.block_bc_flag = block_bc_flag self.material = material self.bc = bc - print 'HEX Blocks:' + print('HEX Blocks:') for m, f in zip(self.block_mat, self.block_flag): - print 'block ', m, 'material flag ', f - print 'Absorbing Boundary Conditions:' + print('block ', m, 'material flag ', f) + print('Absorbing Boundary Conditions:') for m, f in zip(self.block_bc, self.block_bc_flag): - print 'bc ', m, 'bc flag ', f - print 'Topography (free surface)' - print self.topography - print 'Free surface' - print self.free + print('bc ', m, 'bc flag ', f) + print('Topography (free surface)') + print(self.topography) + print('Free surface') + print(self.free) 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 '****************************************' + 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('****************************************') def get_hex_connectivity(self, ind): if self.hex27: @@ -681,7 +682,7 @@ def mat_parameter(self, properties): # format nummaterials file: #material_domain_id #material_id #rho #vp # #vs #Q_kappa #Q_mu #anisotropy_flag flag = properties[1] - print 'number of material:', flag + print('number of material:', flag) if flag > 0: vel = properties[2] if properties[2] is None and type(vel) != str: @@ -740,11 +741,11 @@ def mat_parameter(self, properties): txt = '%1i %3i %s %s \n' % (properties[0], properties[ 1], properties[2], helpstring) # - # print txt + # print(txt) return txt def nummaterial_write(self, nummaterial_name, placeholder=True): - print 'Writing ' + 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) @@ -789,7 +790,7 @@ def nummaterial_write(self, nummaterial_name, placeholder=True): ''' nummaterial.write(txt) nummaterial.close() - print 'Ok' + print('Ok') def create_hexnode_string(self, hexa, hexnode_string=True): nodes = self.get_hex_connectivity(hexa) @@ -798,9 +799,17 @@ def create_hexnode_string(self, hexa, hexnode_string=True): ordered_nodes = [hexa] + list(nodes[:20]) + [nodes[21]] +\ [nodes[25]] + [nodes[24]] + [nodes[26]] +\ [nodes[23]] + [nodes[22]] + [nodes[20]] + #debug + #print('debug: hex27 connectivity original:',nodes) + # connectivity: (4, 7, 8, 3, 1, 6, 5, 2, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27) + #print('debug: hex27 connectivity ordered :',ordered_nodes) + # connectivity: [1, 4, 7, 8, 3, 1, 6, 5, 2, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 22, 26, 25, 27, 24, 23, 21] txt = ' '.join(str(x) for x in ordered_nodes) txt = txt + '\n' else: + #debug + #print('debug: hex8 connectivity original:',nodes) + # connectivity used is: 4, 7, 8, 3, 1, 6, 5, 2 txt = str(hexa) + ' ' + ' '.join(str(x) for x in nodes) txt = txt + '\n' if hexnode_string: @@ -834,30 +843,30 @@ def create_facenode_string(self, hexa, face, normal=None, cknormal=True, map(int, txt.split()) def mesh_write(self, mesh_name): - print 'Writing ' + mesh_name + '..... v2' + print('Writing ' + mesh_name + '..... v2') num_elems = cubit.get_hex_count() meshfile = open(mesh_name, 'w') - print ' total number of elements:', str(num_elems) + print(' total number of elements:', str(num_elems)) meshfile.write(str(num_elems) + '\n') for block, flag in zip(self.block_mat, self.block_flag): hexes = cubit.get_block_hexes(block) - print 'block ', block, ' hexes ', len(hexes) + print('block ', block, ' hexes ', len(hexes)) for hexa in hexes: txt = self.create_hexnode_string(hexa) meshfile.write(txt) meshfile.close() - print 'Ok' + print('Ok') def material_write(self, mat_name): mat = open(mat_name, 'w') - print 'Writing ' + mat_name + '.....' + print('Writing ' + mat_name + '.....') for block, flag in zip(self.block_mat, self.block_flag): - print 'block ', block, 'flag ', flag + print('block ', block, 'flag ', flag) hexes = cubit.get_block_hexes(block) for hexa in hexes: mat.write(('%10i %10i\n') % (hexa, flag)) mat.close() - print 'Ok' + print('Ok') def get_extreme(self, c, cmin, cmax): if not cmin and not cmax: @@ -872,10 +881,10 @@ def get_extreme(self, c, cmin, cmax): def nodescoord_write(self, nodecoord_name): nodecoord = open(nodecoord_name, 'w') - print 'Writing ' + nodecoord_name + '.....' + print('Writing ' + nodecoord_name + '.....') node_list = cubit.parse_cubit_list('node', 'all') num_nodes = len(node_list) - print ' number of nodes:', str(num_nodes) + print(' number of nodes:', str(num_nodes)) nodecoord.write('%10i\n' % num_nodes) # for node in node_list: @@ -886,7 +895,7 @@ def nodescoord_write(self, nodecoord_name): txt = ('%10i %20f %20f %20f\n') % (node, x, y, z) nodecoord.write(txt) nodecoord.close() - print 'Ok' + print('Ok') def free_write(self, freename=None): # free surface @@ -897,17 +906,17 @@ def free_write(self, freename=None): if not freename: freename = self.freename # writes free surface file - print 'Writing ' + freename + '.....' + 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 'free surface (topography) block name:', \ - name, 'id:', block + print('free surface (topography) block name:', \ + name, 'id:', block) quads_all = cubit.get_block_faces(block) - print ' number of faces = ', len(quads_all) + print(' number of faces = ', 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') @@ -920,9 +929,9 @@ def free_write(self, freename=None): freehex.write(txt) elif block == self.free: name = cubit.get_exodus_entity_name('block', block) - print 'free surface block name:', name, 'id:', block + print('free surface block name:', name, 'id:', block) quads_all = cubit.get_block_faces(block) - print ' number of faces = ', len(quads_all) + print(' number of faces = ', 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') @@ -934,7 +943,7 @@ def free_write(self, freename=None): h, f, normal, cknormal=False) freehex.write(txt) freehex.close() - print 'Ok' + print('Ok') cubit.cmd('set info on') cubit.cmd('set echo on') @@ -950,11 +959,11 @@ def check_cmpl_size(self, case='x'): vmintmp = self.zmin if self.size > .3 * (vmaxtmp - vmintmp): - print 'please select the size of cpml less than 30% of the ' +\ - case + ' size of the volume' - print vmaxtmp - vmintmp, .3 * (vmaxtmp - vmintmp) - print 'cmpl set to false, no ' + self.cpmlname +\ - ' file will be created' + print('please select the size of cpml less than 30% of the ' +\ + case + ' size of the volume') + print(vmaxtmp - vmintmp, .3 * (vmaxtmp - vmintmp)) + print('cmpl set to false, no ' + self.cpmlname +\ + ' file will be created') return False, False else: vmin = vmintmp + self.size @@ -1035,10 +1044,10 @@ def abs_write(self, absname=None): if self.cpml: if not absname: absname = self.cpmlname - print 'Writing cpml' + absname + '.....' + print('Writing cpml' + absname + '.....') list_cpml = self.select_cpml() if list_cpml is False: - print 'error writing cpml files' + print('error writing cpml files') return else: abshex_cpml = open(absname, 'w') @@ -1059,57 +1068,57 @@ def abs_write(self, absname=None): 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 + print(' block name:', name, 'id:', block) cknormal = True abshex_local = False # opens file if re.search('xmin', name): - print 'xmin' + print("xmin") abshex_local = open(absname + '_xmin', 'w') normal = (-1, 0, 0) elif re.search('xmax', name): - print "xmax" + print("xmax") abshex_local = open(absname + '_xmax', 'w') normal = (1, 0, 0) elif re.search('ymin', name): - print "ymin" + print("ymin") abshex_local = open(absname + '_ymin', 'w') normal = (0, -1, 0) elif re.search('ymax', name): - print "ymax" + print("ymax") abshex_local = open(absname + '_ymax', 'w') normal = (0, 1, 0) elif re.search('bottom', name): - print "bottom" + print("bottom") abshex_local = open(absname + '_bottom', 'w') normal = (0, 0, -1) elif re.search('abs', name): - print "abs all - experimental, check the output" + print("abs all - experimental, check the output") cknormal = False abshex_local = open(absname, 'w') else: if block == 1003: - print 'xmin' + print("xmin") abshex_local = open(absname + '_xmin', 'w') normal = (-1, 0, 0) elif block == 1004: - print "ymin" + print("ymin") abshex_local = open(absname + '_ymin', 'w') normal = (0, -1, 0) elif block == 1005: - print "xmax" + print("xmax") abshex_local = open(absname + '_xmax', 'w') normal = (1, 0, 0) elif block == 1006: - print "ymax" + print("ymax") abshex_local = open(absname + '_ymax', 'w') normal = (0, 1, 0) elif block == 1002: - print "bottom" + print("bottom") abshex_local = open(absname + '_bottom', 'w') normal = (0, 0, -1) elif block == 1000: - print "custumized" + print("custumized") abshex_local = open(absname, 'w') cknormal = False normal = None @@ -1119,7 +1128,7 @@ def abs_write(self, absname=None): # gets face elements quads_all = cubit.get_block_faces(block) dic_quads_all = dict(zip(quads_all, quads_all)) - print ' number of faces = ', len(quads_all) + print(' number of faces = ', len(quads_all)) abshex_local.write('%10i\n' % len(quads_all)) for h in list_hex: faces = cubit.get_sub_elements('hex', h, 2) @@ -1129,7 +1138,7 @@ def abs_write(self, absname=None): h, f, normal=normal, cknormal=cknormal) abshex_local.write(txt) abshex_local.close() - print 'Ok' + print('Ok') cubit.cmd('set info on') cubit.cmd('set echo on') @@ -1152,13 +1161,13 @@ def surface_write(self, pathdir=None): else: continue # gets face elements - print ' surface block name: ', name, 'id: ', block + print(' surface block name: ', name, 'id: ', block) quads_all = cubit.get_block_faces(block) - print ' face = ', len(quads_all) + print(' face = ', len(quads_all)) if len(quads_all) == 0: continue # writes out surface infos to file - print 'Writing ' + filename + '.....' + print('Writing ' + filename + '.....') surfhex_local = open(filename, 'w') dic_quads_all = dict(zip(quads_all, quads_all)) # writes number of surface elements @@ -1174,17 +1183,17 @@ def surface_write(self, pathdir=None): surfhex_local.write(txt) # closes file surfhex_local.close() - print 'Ok' + print('Ok') def rec_write(self, recname): - print 'Writing ' + 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' + print('Ok') def write(self, path='', netcdf_name=False): cubit.cmd('set info off') @@ -1261,8 +1270,8 @@ def write(self, path='', netcdf_name=False): # if block != self.topography and block != self.free: # label, normal, cknormal = self._get_bc_flag(block) # quads_all = cubit.get_block_faces(block) -# print label -# print ' number of faces = ', len(quads_all) +# print(label) +# print(' number of faces = ', len(quads_all)) # self.netcdf_db.createDimension( # 'num_el_' + label, len(quads_all)) # self.netcdf_db.createVariable( @@ -1285,7 +1294,7 @@ def write(self, path='', netcdf_name=False): def _get_bc_flag(self, block): import re name = cubit.get_exodus_entity_name('block', block) - print ' block name:', name, 'id:', block + print(' block name:', name, 'id:', block) cknormal = True if re.search('xmin', name): label = 'xmin' @@ -1304,7 +1313,7 @@ def _get_bc_flag(self, block): normal = (0, 0, -1) elif re.search('abs', name): label = "all" - print "abs all - experimental, check the output" + print("abs all - experimental, check the output") cknormal = False else: if block == 1003: @@ -1343,16 +1352,16 @@ def _get_bc_flag(self, block): # abs_array.append(self.create_facenode_string( # h, f, normal=normal, cknormal=cknormal, # facenode_string=False)) -# print 'Ok' +# print('Ok') # cubit.cd('set info on') # cubit.cmd('set echo on') # return abs_array # def _netcdf_nodescoord_array(self): -# print 'Writing node coordinates..... netcdf' +# print('Writing node coordinates..... netcdf') # node_list = cubit.parse_cubit_list('node', 'all') # num_nodes = len(node_list) -# print ' number of nodes:', str(num_nodes) +# print(' number of nodes:', str(num_nodes)) # # # coord_array = [] # map_array = [] @@ -1363,7 +1372,7 @@ def _get_bc_flag(self, block): # self.zmin, self.zmax = self.get_extreme(z, self.zmin, self.zmax) # map_array.append([map]) # coord_array.append([x, y, z]) -# print 'Ok' +# print('Ok') # return map_array, coord_array # def _netcdf_free_array(self): @@ -1373,31 +1382,31 @@ def _get_bc_flag(self, block): # cubit.cmd('set journal off') # normal = (0, 0, 1) # # writes free surface file -# print 'Writing free surface..... netcdf' +# print('Writing free surface..... netcdf') # # # # searches block definition with name face_topo # free_array = [] # 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 'free surface (topography) block name:',name,'id:',block +# print('free surface (topography) block name:',name,'id:',block) # quads_all = cubit.get_block_faces(block) -# print ' number of faces = ', len(quads_all) +# print(' number of faces = ', len(quads_all)) # dic_quads_all = dict(zip(quads_all, 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 +# # print(f) # free_array.append(self.create_facenode_string( # h, f, normal, cknormal=True, # facenode_string=False)) # elif block == self.free: # name = cubit.get_exodus_entity_name('block', block) -# print 'free surface block name:', name, 'id:', block +# print('free surface block name:', name, 'id:', block) # quads_all = cubit.get_block_faces(block) -# print ' number of faces = ', len(quads_all) +# print(' number of faces = ', len(quads_all)) # dic_quads_all = dict(zip(quads_all, quads_all)) # list_hex = cubit.parse_cubit_list('hex', 'all') # for h in list_hex: @@ -1407,33 +1416,33 @@ def _get_bc_flag(self, block): # free_array.append(self.create_facenode_string( # h, f, normal, cknormal=False, # facenode_string=False)) -# print 'Ok' +# print('Ok') # cubit.cmd('set info on') # cubit.cmd('set echo on') # return free_array # def _netcdf_material_array(self): -# print 'Writing material...... netcdf' +# print('Writing material...... netcdf') # material_array = [] # for block, flag in zip(self.block_mat, self.block_flag): -# print 'block ', block, 'flag ', flag +# print('block ', block, 'flag ', flag) # hexes = cubit.get_block_hexes(block) # for hexa in hexes: # material_array.append([hexa, flag]) -# print 'Ok' +# print('Ok') # return material_array # def _netcdf_mesh_array(self): -# print 'Writing ' + mesh_name + '..... netcdf' -# print 'total number of elements:', str(cubit.get_hex_count()) +# print('Writing ' + mesh_name + '..... netcdf') +# print('total number of elements:', str(cubit.get_hex_count())) # mesh_array = [] # for block, flag in zip(self.block_mat, self.block_flag): # hexes = cubit.get_block_hexes(block) -# print 'block ', block, ' hexes ', len(hexes) +# print('block ', block, ' hexes ', len(hexes)) # for hexa in hexes: # mesh_array.append(create_hexnode_string( # hexa, hexnode_string=False)) -# print 'Ok' +# print('Ok') # return mesh_array # def _write_ascii(self, path=''): @@ -1463,11 +1472,11 @@ def export2SPECFEM3D(path_exporting_mesh_SPECFEM3D='.', hex27=False, cpml=False, cpml_size=False, top_absorbing=False): sem_mesh = mesh(hex27, cpml, cpml_size, top_absorbing) # sem_mesh.block_definition() - # print sem_mesh.block_mat - # print sem_mesh.block_flag + # print(sem_mesh.block_mat) + # print(sem_mesh.block_flag) # sem_mesh.write(path=path_exporting_mesh_SPECFEM3D) - print 'END SPECFEM3D exporting process......' + print('END SPECFEM3D exporting process......') if cpml: cmd = 'save as "cpml.cub" overwrite' cubit.cmd(cmd) diff --git a/CUBIT_GEOCUBIT/geocubitlib/dev/cmpl.py b/CUBIT_GEOCUBIT/geocubitlib/dev/cmpl.py index 2ac583ae5..b10ea7cdd 100755 --- a/CUBIT_GEOCUBIT/geocubitlib/dev/cmpl.py +++ b/CUBIT_GEOCUBIT/geocubitlib/dev/cmpl.py @@ -23,7 +23,8 @@ # 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. # # # ############################################################################# -# +from __future__ import print_function + try: import start as start cubit = start.start_cubit() @@ -31,7 +32,7 @@ try: import cubit except: - print 'error importing cubit, check if cubit is installed' + print('error importing cubit, check if cubit is installed') pass @@ -119,7 +120,7 @@ def mesh_cpml(list_vol,remesh=True,refinement=None,top_surf=None,size=None): for refdepth in refinement: cubit.cmd('refine surf '+top_surf+' numsplit 1 bias 1 depth '+str(refdepth)) except: - print 'DEBUG: error in refinement cpml' + print('DEBUG: error in refinement cpml') xmin=xmin-size xmax=xmax+size ymin=ymin-size @@ -145,7 +146,7 @@ def collecting_cpml(ip,size=None,cpuxmin=0,cpuxmax=1,cpuymin=0,cpuymax=1,cpux=1, from utilities import list2str # if not size: - print 'cpml size must be specified' + print('cpml size must be specified') return @@ -187,7 +188,7 @@ def collecting_cpml(ip,size=None,cpuxmin=0,cpuxmax=1,cpuymin=0,cpuymax=1,cpux=1, # # - #print boundary_dict + #print(boundary_dict) block_list=cubit.get_block_id_list() for block in block_list: ty=cubit.get_block_element_type(block) diff --git a/CUBIT_GEOCUBIT/geocubitlib/exportlib.py b/CUBIT_GEOCUBIT/geocubitlib/exportlib.py index 11e1f3187..7e301ad0b 100755 --- a/CUBIT_GEOCUBIT/geocubitlib/exportlib.py +++ b/CUBIT_GEOCUBIT/geocubitlib/exportlib.py @@ -23,7 +23,8 @@ # 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. # # # ############################################################################# -# +from __future__ import print_function + try: import start as start cubit = start.start_cubit() @@ -31,12 +32,12 @@ try: import cubit except: - print "error importing cubit, check if cubit is installed" + print("error importing cubit, check if cubit is installed") pass import glob from utilities import get_cubit_version -print 'version 2.2' +print('version 2.2') class VersionException(Exception): @@ -107,11 +108,11 @@ def collect_new(cpuxmin=0, cpuxmax=1, cpuymin=0, cpuymax=1, cpux=1, cpuy=1, # cubit.cmd('set error off') version_cubit = get_cubit_version() - print '##########################################################' - print '#' - print '# collecting mesh' - print '#' - print '##########################################################' + print('##########################################################') + print('#') + print('# collecting mesh') + print('#') + print('##########################################################') if decimate: if version_cubit >= 14.0: @@ -205,11 +206,11 @@ def collect_new(cpuxmin=0, cpuxmax=1, cpuymin=0, cpuymax=1, cpux=1, cpuy=1, cubit.cmd('set info echo journ on') # # - print 'end meshing' + print('end meshing') # # if qlog: - print '\n\nQUALITY CHECK.... ***************\n\n' + print('\n\nQUALITY CHECK.... ***************\n\n') import quality_log tq = open(outdir2 + outfilename + '.quality', 'w') max_skewness, min_length = quality_log.quality_log(tq) @@ -228,7 +229,7 @@ def collect_new(cpuxmin=0, cpuxmax=1, cpuymin=0, cpuymax=1, cpux=1, cpuy=1, ' '.join(str(x) for x in vol_blocks) + " feature_angle 135.0") command = "save as '" + outdir2 + outfilename + ".cub' overwrite" - print command + print(command) cubit.cmd(command) @@ -238,14 +239,14 @@ def e2SEM(files=False, listblock=None, listflag=None, outdir='.', if files: filenames = glob.glob(files) for f in filenames: - print f + print(f) extension = f.split('.')[-1] if extension == 'cub': cubit.cmd('open "' + f + '"') elif extension == 'e': cubit.cmd('import mesh "' + f + '" no_geom') else: - print extension + print(extension) if listblock and listflag: pass else: @@ -274,14 +275,14 @@ def collecting_block(store_group_name, iproc=0, xmin=[0], xmax=[0], ymin=[0], block_list.sort() block_hex = [x for x in block_list if x <= 1000] block_side = [x for x in block_list if x > 1000] - #print 'collecting block list = ',block_list - #print ' hex list = ',block_hex - #print ' side list = ',block_side - print 'collecting block: iproc = ',iproc,xmin,xmax,ymin,ymax - print ' block_hex = ',block_hex,'block_side = ',block_side + #print('collecting block list = ',block_list) + #print(' hex list = ',block_hex) + #print(' side list = ',block_side) + print('collecting block: iproc = ',iproc,xmin,xmax,ymin,ymax) + print(' block_hex = ',block_hex,'block_side = ',block_side) ## hexahedrals for ib, block in enumerate(block_hex): - #print ' hex ib = ',ib,'block = ',block + #print(' hex ib = ',ib,'block = ',block) if index_block == 0: cubit.cmd("group 'vol" + str(block) + "' add Hex in block " + str(block)) @@ -293,7 +294,7 @@ def collecting_block(store_group_name, iproc=0, xmin=[0], xmax=[0], ymin=[0], cubit.cmd("del block " + str(block)) ## faces for ib, side in enumerate(block_side): - #print ' faces ib = ',ib,'side = ',side + #print(' faces ib = ',ib,'side = ',side) if side == 1004: if iproc in ymin: cubit.cmd("group 'ymin' add face in block " + str(side)) @@ -323,7 +324,7 @@ def collecting_block(store_group_name, iproc=0, xmin=[0], xmax=[0], ymin=[0], ## check if faces in group lateral ilateral = cubit.get_id_from_name('lateral') lateral_nodes = cubit.get_group_nodes(ilateral) - print ' lateral nodes: ', len(lateral_nodes) + print(' lateral nodes: ', len(lateral_nodes)) return store_group_name @@ -333,12 +334,12 @@ def check_lateral_nodes(name_group='lateral'): ilateral_nodes = cubit.get_id_from_name('lateral_nodes') lateral_nodes = cubit.get_group_nodes(ilateral_nodes) cubit.cmd('del group ' + str(ilateral_nodes)) - print name_group, ' nodes ', len(lateral_nodes) + print(name_group, ' nodes ', len(lateral_nodes)) return lateral_nodes def prepare_equivalence_new(name_group='lateral'): - print 'equivalence group ',name_group + print('equivalence group ',name_group) length = {} cmd = "group 'tmpn' add edge in face in group " + name_group cubit.cmd(cmd) @@ -350,7 +351,7 @@ def prepare_equivalence_new(name_group='lateral'): length[e] = lengthmin * 0.5 cubit.cmd('delete group ' + str(ge)) - #print ' equivalence edge lengths ',length + #print(' equivalence edge lengths ',length) if len(length) > 0: minvalue = min(length.values()) maxvalue = max(length.values()) @@ -358,7 +359,7 @@ def prepare_equivalence_new(name_group='lateral'): minvalue = 100. maxvalue = 100. # - print ' min length: ', minvalue, 'max length: ', maxvalue + print(' min length: ', minvalue, 'max length: ', maxvalue) if minvalue != 0: nbin = int((maxvalue / minvalue)) + 1 factor = (maxvalue - minvalue) / nbin @@ -378,26 +379,26 @@ def prepare_equivalence_new(name_group='lateral'): for k in range(0, len(inv_length.keys()) - 1): inv_length[ks[k]] = inv_length[ks[k]] + inv_length[ks[k + 1]] - print ' edge lengths ',inv_length.keys(), factor, minvalue + print(' edge lengths ',inv_length.keys(), factor, minvalue) return factor, minvalue, maxvalue, inv_length def merging_node_new(tol, clean=True, graphic_debug=False): empty = False - print 'tolerance ', tol + print('tolerance ', tol) cubit.cmd("topology check coincident node node in \ group coincident_lateral_nodes tolerance " + str(tol) + " highlight brief result \ group 'merging_lateral_nodes'") group_exist = cubit.get_id_from_name("merging_lateral_nodes") if not group_exist: - print 'no nodes in this tolerance range' + print('no nodes in this tolerance range') else: merging_nodes = cubit.get_group_nodes(group_exist) if graphic_debug: cubit.cmd('draw group lateral') cubit.cmd('high group merging_lateral_nodes') - print 'merging ', len(merging_nodes), ' nodes.....' + print('merging ', len(merging_nodes), ' nodes.....') cubit.cmd("equivalence node in merging_lateral_nodes \ tolerance " + str(tol * 2)) if clean: @@ -406,7 +407,7 @@ def merging_node_new(tol, clean=True, graphic_debug=False): cubit.cmd("delete Group merging_lateral_nodes") ic_nodes = cubit.get_id_from_name('coincident_lateral_nodes') c_nodes = cubit.get_group_nodes(ic_nodes) - print len(c_nodes) + print(len(c_nodes)) if len(c_nodes) == 0: empty = True if graphic_debug: @@ -436,8 +437,8 @@ def graphic_merging(tol, step_tol=None, maxtol=None): isempty = merging_node_new(tol, clean=True, graphic_debug=True) tol = tol + step_tol if tol > maxtol: - print 'tolerance greater than the max length of the edges, \ - please check the mesh' + print('tolerance greater than the max length of the edges, \ + please check the mesh') def collecting_merging_new(cpuxmin=0, cpuxmax=0, cpuymin=0, cpuymax=0, cpux=1, @@ -453,27 +454,27 @@ def collecting_merging_new(cpuxmin=0, cpuxmax=0, cpuymin=0, cpuymax=0, cpux=1, pass # number_of_chunks = cpux * cpuy - print 'number of chunks: ', number_of_chunks + print('number of chunks: ', number_of_chunks) xmin, xmax, ymin, ymax, listfull = map_boundary(cpuxmin, cpuxmax, cpuymin, cpuymax, cpux, cpuy) - print 'list of processes with boundary:' - print ' xmin: ', xmin - print ' xmax: ', xmax - print ' ymin: ', ymin - print ' ymax: ', ymax - print ' full list: ', listfull + print('list of processes with boundary:') + print(' xmin: ', xmin) + print(' xmax: ', xmax) + print(' ymin: ', ymin) + print(' ymax: ', ymax) + print(' full list: ', listfull) if 1 < number_of_chunks < max(listfull): raise MergingError('error mapping the chunks') # if cubfiles: nf, listiproc, filenames, cubflag = importing_cubfiles(cubfiles) - print nf, listiproc, filenames, cubflag, listfull + print(nf, listiproc, filenames, cubflag, listfull) else: nf = 0 filenames = [] iproc = 0 - print 'cubfiles : ',cubfiles - print 'number of files: ',nf + print('cubfiles : ',cubfiles) + print('number of files: ',nf) # index_block = -1 @@ -488,10 +489,10 @@ def collecting_merging_new(cpuxmin=0, cpuxmax=0, cpuymin=0, cpuymax=0, cpux=1, if nf > 0: for iproc, filename in zip(listiproc, filenames): - print iproc, filename, iproc in listfull + print(iproc, filename, iproc in listfull) try: if iproc in listfull: - print filename + print(filename) index_block = index_block + 1 if cubflag: cubit.cmd('import cubit "' + filename + '"') @@ -499,15 +500,15 @@ def collecting_merging_new(cpuxmin=0, cpuxmax=0, cpuymin=0, cpuymax=0, cpux=1, cubit.cmd('import mesh "' + filename + '" block all no_geom') except: cubit.cmd('import mesh "' + filename + '" block all no_geom') - # print iproc,xmin,xmax,ymin,ymax,iproc in xmin,iproc in xmax,iproc in ymin,iproc in ymax - print "iproc ",iproc," mesh imported" + # print(iproc,xmin,xmax,ymin,ymax,iproc in xmin,iproc in xmax,iproc in ymin,iproc in ymax) + print("iproc ",iproc," mesh imported") store_tmp = collecting_block(store_group_name, iproc, xmin, xmax, ymin, ymax, index_block) - print 'store tmp ',store_tmp + print('store tmp ',store_tmp) if len(store_tmp) != 0: store_group_name = store_tmp # lateral_nodes = check_lateral_nodes() - print 'store group name ',store_group_name + print('store group name ',store_group_name) cubit.cmd('save as "tmp_nomerging.cub" overwrite ') @@ -555,10 +556,10 @@ def collecting_merging_new(cpuxmin=0, cpuxmax=0, cpuymin=0, cpuymax=0, cpux=1, for ig, g in enumerate(store_group_name): cubit.cmd('block ' + str(ig + 1) + ' hex in group ' + g) cubit.cmd('block ' + str(ig + 1) + ' name "vol' + str(ig + 1) + '"') - print 'block ' + str(ig + 1) + ' hex in group ' + g + print('block ' + str(ig + 1) + ' hex in group ' + g) for ig, g in enumerate(side_name): cubit.cmd('block ' + side_val[ig] + ' face in group ' + g) - print 'block ' + side_val[ig] + ' face in group ' + g + print('block ' + side_val[ig] + ' face in group ' + g) cubit.cmd('block ' + side_val[ig] + ' name "' + side_block_name[ig] + '"') cubit.cmd('del group all') @@ -660,10 +661,10 @@ def refine_closecurve(block=1001, closed_filenames=None, acis=True): if not isinstance(closed_filenames, list): closed_filenames = [closed_filenames] for f in closed_filenames: - print f + print(f) if acis: curves = curves + load_curves(f) - print curves + print(curves) blist = list(cubit.get_block_id_list()) try: blist.remove(1001) @@ -838,7 +839,7 @@ def collecting_merging(cpuxmin=0, cpuxmax=1, cpuymin=0, cpuymax=1, cpux=1, cpuy, cpuxmin, cpuxmax, cpuymin, cpuymax) # # - # print boundary_dict + # print(boundary_dict) block_list = cubit.get_block_id_list() for block in block_list: ty = cubit.get_block_element_type(block) @@ -846,7 +847,7 @@ def collecting_merging(cpuxmin=0, cpuxmax=1, cpuymin=0, cpuymax=1, cpux=1, cubit.cmd('block ' + str(block) + ' name "vol' + str(block) + '"') # # - print 'chbound', ckbound_method1, ckbound_method2 + print('chbound', ckbound_method1, ckbound_method2) if ckbound_method1 and not ckbound_method2 and len(filenames) != 1: # use the equivalence method for groups @@ -866,7 +867,7 @@ def collecting_merging(cpuxmin=0, cpuxmax=1, cpuymin=0, cpuymax=1, cpux=1, for iy in range(cpuymin, cpuymax): ind = ind + 1 iproc = iy * cpux + ix - print '******************* ', iproc, ind, '/', len(listfull) + print('******************* ', iproc, ind, '/', len(listfull)) # # ileft | iproc # -------------------- @@ -890,7 +891,7 @@ def collecting_merging(cpuxmin=0, cpuxmax=1, cpuymin=0, cpuymax=1, cpux=1, idown = iproc idiag = ileft # - print iproc, ileft, idiag, idown + print(iproc, ileft, idiag, idown) if iproc != idown: nup = boundary_dict[iproc]['nodes_surf_ymin'] ndow = boundary_dict[idown]['nodes_surf_ymax'] @@ -952,8 +953,8 @@ def collecting_merging(cpuxmin=0, cpuxmax=1, cpuymin=0, cpuymax=1, cpux=1, group_id_1 = cubit.get_id_from_name("negativejac") n1 = cubit.get_group_nodes(group_id_1) if len(n1) != 0: - print 'error, negative jacobian after the equivalence node command, \ - use --merge2 instead of --equivalence/--merge/--merge1' + print('error, negative jacobian after the equivalence node command, \ + use --merge2 instead of --equivalence/--merge/--merge1') elif ckbound_method2 and not ckbound_method1 and len(filenames) != 1: if isinstance(merge_tolerance, list): tol = merge_tolerance[0] @@ -966,7 +967,7 @@ def collecting_merging(cpuxmin=0, cpuxmax=1, cpuymin=0, cpuymax=1, cpux=1, for ix in range(cpuxmin, cpuxmax): for iy in range(cpuymin, cpuymax): iproc = iy * cpux + ix - print '******************* ', iproc + print('******************* ', iproc) # # ileft | iproc # -------------------- @@ -1059,8 +1060,8 @@ def collecting_merging(cpuxmin=0, cpuxmax=1, cpuymin=0, cpuymax=1, cpux=1, if group_id_1 != 0: n1 = cubit.get_group_nodes(group_id_1) if len(n1) != 0: - print 'error, coincident nodes after the equivalence \ - node command, check the tolerance' + print('error, coincident nodes after the equivalence \ + node command, check the tolerance') import sys sys.exit() cmd = 'group "negativejac" add quality hex all Jacobian high' @@ -1068,8 +1069,8 @@ def collecting_merging(cpuxmin=0, cpuxmax=1, cpuymin=0, cpuymax=1, cpux=1, group_id_1 = cubit.get_id_from_name("negativejac") n1 = cubit.get_group_nodes(group_id_1) if len(n1) != 0: - print 'error, negative jacobian after the equivalence node command, \ - check the mesh' + print('error, negative jacobian after the equivalence node command, \ + check the mesh') elif ckbound_method1 and ckbound_method2 and len(filenames) != 1: block_list = cubit.get_block_id_list() i = -1 @@ -1090,11 +1091,11 @@ def collecting_merging(cpuxmin=0, cpuxmax=1, cpuymin=0, cpuymax=1, cpux=1, str(block) + ' tolerance ' + str(tol) + \ ' nodraw brief result group "b' + str(block) + '"' cubit.cmd(cmd) - print cmd + print(cmd) cmd = 'equivalence node in group b' + \ str(block) + ' tolerance ' + str(tol) cubit.cmd(cmd) - print cmd + print(cmd) if isinstance(merge_tolerance, list): tol = max(merge_tolerance) elif merge_tolerance: @@ -1110,8 +1111,8 @@ def collecting_merging(cpuxmin=0, cpuxmax=1, cpuymin=0, cpuymax=1, cpux=1, if group_id_1 != 0: n1 = cubit.get_group_nodes(group_id_1) if len(n1) != 0: - print 'error, coincident nodes after the equivalence node \ - command, check the tolerance' + print('error, coincident nodes after the equivalence node \ + command, check the tolerance') import sys sys.exit() cmd = 'group "negativejac" add quality hex all Jacobian high' @@ -1119,8 +1120,8 @@ def collecting_merging(cpuxmin=0, cpuxmax=1, cpuymin=0, cpuymax=1, cpux=1, group_id_1 = cubit.get_id_from_name("negativejac") n1 = cubit.get_group_nodes(group_id_1) if len(n1) != 0: - print 'error, negative jacobian after the equivalence node command, \ - use --merge instead of --equivalence' + print('error, negative jacobian after the equivalence node command, \ + use --merge instead of --equivalence') def collect_old(cpuxmin=0, cpuxmax=1, cpuymin=0, cpuymax=1, cpux=1, cpuy=1, @@ -1189,11 +1190,11 @@ def collect_old(cpuxmin=0, cpuxmax=1, cpuymin=0, cpuymax=1, cpux=1, cpuy=1, command = "save as '" + outdir2 + outfilename + ".cub' overwrite" cubit.cmd(command) # - print 'end meshing' + print('end meshing') # # if qlog: - print '\n\nQUALITY CHECK.... ***************\n\n' + print('\n\nQUALITY CHECK.... ***************\n\n') import quality_log tq = open(outdir2 + outfilename + '.quality', 'w') max_skewness, min_length = quality_log.quality_log(tq) @@ -1213,14 +1214,14 @@ def e2SEM_old(files=False, listblock=None, listflag=None, outdir='.', if files: filenames = glob.glob(files) for f in filenames: - print f + print(f) extension = f.split('.')[-1] if extension == 'cub': cubit.cmd('open "' + f + '"') elif extension == 'e': cubit.cmd('import mesh "' + f + '" no_geom') else: - print extension + print(extension) if listblock and listflag: pass else: @@ -1263,7 +1264,7 @@ def prepare_equivalence(nodes1, nodes2): cubit.cmd('delete group ' + str(ge)) minvalue = min(length.values()) maxvalue = max(length.values()) - print 'min lentgh: ', minvalue, 'max lentgh: ', maxvalue + print('min lentgh: ', minvalue, 'max lentgh: ', maxvalue) nbin = int((maxvalue / minvalue) / 2.) + 1 factor = (maxvalue - minvalue) / nbin dic_new = {} @@ -1273,7 +1274,7 @@ def prepare_equivalence(nodes1, nodes2): else: dic_new[k] = 0. inv_length = invert_dict(dic_new) - print inv_length.keys(), factor, minvalue + print(inv_length.keys(), factor, minvalue) ks = inv_length.keys() ks.sort() for k in range(0, len(inv_length.keys()) - 1): @@ -1299,9 +1300,9 @@ def merge_node_ck(n1, n2): for x in inv_length[k]) +\ ' tolerance ' + str(k * factor + minvalue / 3.) cubit.cmd(cmd) - print 'equivalence ' + str(len(inv_length[k])) +\ + print('equivalence ' + str(len(inv_length[k])) +\ ' couples of nodes - tolerance ' + \ - str(k * factor + minvalue / 3.) + str(k * factor + minvalue / 3.)) cubit.cmd('group "checkmerge" add node ' + ' '.join(str(n) for n in n1) + @@ -1309,11 +1310,11 @@ def merge_node_ck(n1, n2): idg = cubit.get_id_from_name('checkmerge') remainnodes = cubit.get_group_nodes(idg) - print 'from ' + str(len(n1) + len(n2)) + ' nodes -> ' + \ - str(len(remainnodes)) + ' nodes' + print('from ' + str(len(n1) + len(n2)) + ' nodes -> ' + \ + str(len(remainnodes)) + ' nodes') if len(n1) != len(remainnodes): - print 'equivalence ' + str(len(remainnodes)) + \ - ' couples of nodes - tolerance ' + str(minvalue / 3.) + print('equivalence ' + str(len(remainnodes)) + \ + ' couples of nodes - tolerance ' + str(minvalue / 3.)) cubit.cmd('set info on') cubit.cmd('set echo on') cubit.cmd('set journal on') @@ -1328,7 +1329,7 @@ def merge_node_ck(n1, n2): dimension 3 block all overwrite') cubit.cmd('save as "error_merging.cub" \ dimension 3 block all overwrite') - print 'error merging ' + print('error merging ') if False: import sys sys.exit(2) @@ -1356,9 +1357,9 @@ def merge_node(n1, n2): ' tolerance ' + str(k * factor + minvalue / 3.) cubit.cmd(cmd) - print 'equivalence ' + str(len(inv_length[k])) + \ + print('equivalence ' + str(len(inv_length[k])) + \ ' couples of nodes - tolerance ' + \ - str(k * factor + minvalue / 3.) + str(k * factor + minvalue / 3.)) cubit.cmd('set info on') cubit.cmd('set echo on') @@ -1392,13 +1393,13 @@ def prepare_equivalence_4(nodes1, nodes2, nodes3, nodes4): maxvalue = max(length.values()) except: try: - print nodes - print 'edges ', e1 + print(nodes) + print('edges ', e1) except: pass minvalue = 10. maxvalue = 2000. - print 'min lentgh: ', minvalue, 'max lentgh: ', maxvalue + print('min lentgh: ', minvalue, 'max lentgh: ', maxvalue) nbin = int((maxvalue / minvalue) / 2.) + 1 factor = (maxvalue - minvalue) / nbin dic_new = {} @@ -1408,7 +1409,7 @@ def prepare_equivalence_4(nodes1, nodes2, nodes3, nodes4): else: dic_new[k] = 0. inv_length = invert_dict(dic_new) - print inv_length.keys(), factor, minvalue + print(inv_length.keys(), factor, minvalue) ks = inv_length.keys() ks.sort() for k in range(0, len(inv_length.keys()) - 1): @@ -1430,7 +1431,7 @@ def get_z(node): def merge_node_4(n1, n2, n3, n4, newmethod=True): if newmethod: - print "merge node 4 side" + print("merge node 4 side") n1o = ording_z(n1) n2o = ording_z(n2) n3o = ording_z(n3) @@ -1455,21 +1456,21 @@ def merge_node_4(n1, n2, n3, n4, newmethod=True): ' '.join(' '.join(str(n) for n in x)) + \ ' tolerance ' + str(k * factor + minvalue / 3.) except: - print k, "***************************************** s" - print inv_length[k] + print(k, "***************************************** s") + print(inv_length[k]) cubit.cmd(cmd) - print 'equivalence ' + str(len(inv_length[k])) +\ + print('equivalence ' + str(len(inv_length[k])) +\ ' couples of nodes - tolerance ' + \ - str(k * factor + minvalue / 3.) + str(k * factor + minvalue / 3.)) if len(inv_length[k]) == 1: cmd = 'equivalence node ' + \ ' '.join(' '.join(str(n) for n in inv_length[k])) + \ ' tolerance ' + str(k * factor + minvalue / 3.) cubit.cmd(cmd) - print 'equivalence ' + str(len(inv_length[k])) + \ + print('equivalence ' + str(len(inv_length[k])) + \ ' couples of nodes - tolerance ' + \ - str(k * factor + minvalue / 3.) + str(k * factor + minvalue / 3.)) ########################################################################################## # diff --git a/CUBIT_GEOCUBIT/geocubitlib/hex_metric.py b/CUBIT_GEOCUBIT/geocubitlib/hex_metric.py index 896577947..aa6147afd 100755 --- a/CUBIT_GEOCUBIT/geocubitlib/hex_metric.py +++ b/CUBIT_GEOCUBIT/geocubitlib/hex_metric.py @@ -22,11 +22,13 @@ # 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. # # # ############################################################################# +from __future__ import print_function + # mesh=SEM_metric_3D # mesh.check_metric() -# print mesh -# -# +# print(mesh) + + import math try: @@ -36,7 +38,7 @@ try: import cubit except: - print 'error importing cubit, check if cubit is installed' + print('error importing cubit, check if cubit is installed') pass @@ -66,22 +68,22 @@ def __repr__(self): if self.max_skewness is not None: self.skew_hystogram = self.hyst( 0, self.max_skewness, self.skew_hyst) - print '-' * 70 - print 'SKEWNESS' - print + print('-' * 70) + print('SKEWNESS') + print('') if len(self.hex_max_skewness) <= 30: - print 'max = ', self.max_skewness, ' in hexes ', \ - self.hex_max_skewness - print '(angle -> minimun =', self.min_angle, \ - ' maximun =', self.max_angle, ')' + print('max = ', self.max_skewness, ' in hexes ', \ + self.hex_max_skewness) + print('(angle -> minimun =', self.min_angle, \ + ' maximun =', self.max_angle, ')') else: - print 'max = ', self.max_skewness, ' in ', \ - len(self.hex_max_skewness), ' hexes ' - print '(angle -> minimun =', self.min_angle, \ - ' maximun =', self.max_angle, ')' - print - print 'skew hystogram' - print + print('max = ', self.max_skewness, ' in ', \ + len(self.hex_max_skewness), ' hexes ') + print('(angle -> minimun =', self.min_angle, \ + ' maximun =', self.max_angle, ')') + print('') + print('skew hystogram') + print('') tot = 0 for i in self.skew_hystogram.values(): tot = tot + len(i) @@ -92,13 +94,13 @@ def __repr__(self): if self.skew_hystogram.has_key(i): if (i + 1) * factor <= 1: nshi = len(self.skew_hystogram[i]) - print i, ' [', i * factor, '->', (i + 1) * factor, \ + print(i, ' [', i * factor, '->', (i + 1) * factor, \ '[ : ', nshi, '/', tot, \ - ' hexes (', nshi / float(tot) * 100., '%)' + ' hexes (', nshi / float(tot) * 100., '%)') else: if (i + 1) * factor <= 1: - print i, ' [', i * factor, '->', (i + 1) * factor, \ - '[ : ', 0, '/', tot, ' hexes (0%)' + print(i, ' [', i * factor, '->', (i + 1) * factor, \ + '[ : ', 0, '/', tot, ' hexes (0%)') print ############################################### if self.min_edge_length is not None: @@ -106,26 +108,26 @@ def __repr__(self): self.min_edge_length, self.max_edge_length, self.edgemin_hyst) self.edgemax_hystogram = self.hyst( self.min_edge_length, self.max_edge_length, self.edgemax_hyst) - print '-' * 70 - print 'edge length' - print + print('-' * 70) + print('edge length') + print('') if len(self.hex_min_edge_length) <= 30: - print 'minimum edge length: ', self.min_edge_length, \ - ' in hexes ', self.hex_min_edge_length + print('minimum edge length: ', self.min_edge_length, \ + ' in hexes ', self.hex_min_edge_length) else: - print 'minimum edge length: ', self.min_edge_length, ' in ', \ - len(self.hex_min_edge_length), ' hexes.' + print('minimum edge length: ', self.min_edge_length, ' in ', \ + len(self.hex_min_edge_length), ' hexes.') if len(self.hex_max_edge_length) <= 30: - print 'maximum edge length: ', self.max_edge_length, \ - ' in hexes ', self.hex_max_edge_length + print('maximum edge length: ', self.max_edge_length, \ + ' in hexes ', self.hex_max_edge_length) else: - print 'maximum edge length: ', self.max_edge_length, \ - ' in ', len(self.hex_max_edge_length), ' hexes.' - print - print 'edge length hystogram' - print + print('maximum edge length: ', self.max_edge_length, \ + ' in ', len(self.hex_max_edge_length), ' hexes.') + print('') + print('edge length hystogram') + print('') factor = (self.max_edge_length - self.min_edge_length) / self.nbin - print 'minimum edge length' + print('minimum edge length') tot = 0 for i in self.edgemin_hystogram.values(): tot = tot + len(i) @@ -134,16 +136,16 @@ def __repr__(self): for i in range(0, self.nbin + 1): if self.edgemin_hystogram.has_key(i): nh = len(self.edgemin_hystogram[i]) - print i, ' [', i * factor + self.min_edge_length, '->', \ - (i + 1) * factor + self.min_edge_length, \ - '[ : ', nh, '/', tot, \ - ' hexes (', nh / float(tot) * 100., '%)' + print(i, ' [', i * factor + self.min_edge_length, '->', \ + (i + 1) * factor + self.min_edge_length, \ + '[ : ', nh, '/', tot, \ + ' hexes (', nh / float(tot) * 100., '%)') else: - print i, ' [', i * factor + self.min_edge_length, '->', \ - (i + 1) * factor + self.min_edge_length, \ - '[ : ', 0, '/', tot, ' hexes (0%)' - print - print 'maximum edge length' + print(i, ' [', i * factor + self.min_edge_length, '->', \ + (i + 1) * factor + self.min_edge_length, \ + '[ : ', 0, '/', tot, ' hexes (0%)') + print('') + print('maximum edge length') tot = 0 for i in self.edgemax_hystogram.values(): tot = tot + len(i) @@ -152,21 +154,21 @@ def __repr__(self): for i in range(0, self.nbin + 1): if self.edgemax_hystogram.has_key(i): nh = len(self.edgemax_hystogram[i]) - print i, ' [', i * factor + self.min_edge_length, \ - '->', (i + 1) * factor + self.min_edge_length, \ - '[ : ', len(self.edgemax_hystogram[i]), '/', tot, \ - ' hexes (', nh / float(tot) * 100., '%)' + print(i, ' [', i * factor + self.min_edge_length, \ + '->', (i + 1) * factor + self.min_edge_length, \ + '[ : ', len(self.edgemax_hystogram[i]), '/', tot, \ + ' hexes (', nh / float(tot) * 100., '%)') else: - print i, ' [', i * factor + self.min_edge_length, \ - '->', (i + 1) * factor + self.min_edge_length, \ - '[ : ', 0, '/', tot, ' hexes (0%)' + print(i, ' [', i * factor + self.min_edge_length, \ + '->', (i + 1) * factor + self.min_edge_length, \ + '[ : ', 0, '/', tot, ' hexes (0%)') if self.dt is not None: - print '-' * 70 - print - print 'STABILITY' - print - print 'time step < ', self.dt, 's, for velocity = ', self.velocity - print + print('-' * 70) + print('') + print('STABILITY') + print('') + print('time step < ', self.dt, 's, for velocity = ', self.velocity) + print('') try: return str(len(self.list_hex)) + ' hexes checked' except: @@ -194,7 +196,7 @@ def hex_metric(self, h): faces = [[0, 1, 2, 3], [4, 5, 6, 7], [1, 5, 6, 2], [0, 4, 7, 3], [2, 6, 7, 3], [1, 5, 4, 0]] else: - print 'bad definition of nodes' + print('bad definition of nodes') return None, None, None x = [] y = [] @@ -218,7 +220,7 @@ def hex_metric(self, h): norm2 = math.sqrt(vx2 * vx2 + vy2 * vy2 + vz2 * vz2) # if norm1 == 0 or norm2 == 0: - print 'degenerated mesh, 0 length edge' + print('degenerated mesh, 0 length edge') import sys sys.exit() angle = math.acos( @@ -284,7 +286,7 @@ def pick_hex(self, list_hex=None, volume=None): cubit.silent_cmd(command) elif list_hex is None and self.list_hex is None: self.list_hex = cubit.parse_cubit_list('hex', 'all') - print 'list_hex: ', len(self.list_hex), ' hexes' + print('list_hex: ', len(self.list_hex), ' hexes') # # ########################################################################## @@ -380,10 +382,10 @@ def check_valence(self): # vmax = max(inv_lookup.keys()) nmax = inv_lookup[max(inv_lookup.keys())] - print 'max hex valence ', vmax, ' at node ', nmax - print '_____' + print('max hex valence ', vmax, ' at node ', nmax) + print('_____') for v in inv_lookup.keys(): - print ('valence %2i - %9i hexes ' % (v, len(inv_lookup[v]))) + print(('valence %2i - %9i hexes ' % (v, len(inv_lookup[v])))) self.valence_summary = inv_lookup self.max_valence = vmax self.node_with_max_valence = nmax @@ -586,7 +588,7 @@ def __init__(self, list_hex=None, volume=None): # def __repr__(self): - print 'check mesh stability' + print('check mesh stability') # def check_simulation_parameter(self, list_hex=None, volume=None, @@ -603,8 +605,8 @@ def check_simulation_parameter(self, list_hex=None, volume=None, # for ind, h in enumerate(self.list_hex): if ind % 10000 == 0: - print 'hex checked: ' + \ - str(int(float(ind) / len(self.list_hex) * 100)) + '%' + print('hex checked: ' + \ + str(int(float(ind) / len(self.list_hex) * 100)) + '%') dt_tmp, pmax_tmp, _, _ = self.hex_simulation_parameter( h, vp_static=vp_static, vs_static=vs_static) timestep_hyst[h] = dt_tmp @@ -622,9 +624,9 @@ def check_simulation_parameter(self, list_hex=None, volume=None, def read_tomo(self, tomofile=None): if tomofile: - print 'reading tomography file ', tomofile + print('reading tomography file ', tomofile) # xtomo,ytomo,ztomo,vp,vs,rho=numpy.loadtxt(tomofile,skiprows=4) - print 'tomography file loaded' + print('tomography file loaded') tf = open(tomofile, 'r') orig_x, orig_y, orig_z, end_x, end_y, end_z = map( float, tf.readline().split()) @@ -648,7 +650,7 @@ def read_tomo(self, tomofile=None): vp.append(v1) vs.append(v2) tf.close() - print 'tomography file loaded' + print('tomography file loaded') # self.orig_x, self.orig_y, self.orig_z, self.end_x, self.end_y, \ self.end_z = orig_x, orig_y, orig_z, end_x, end_y, end_z @@ -660,7 +662,7 @@ def read_tomo(self, tomofile=None): self.xtomo, self.ytomo, self.ztomo, self.vp, self.vs = xtomo, \ ytomo, ztomo, vp, vs else: - print 'no tomofile!!!!' + print('no tomofile!!!!') # # # diff --git a/CUBIT_GEOCUBIT/geocubitlib/local_volume.py b/CUBIT_GEOCUBIT/geocubitlib/local_volume.py index 966f9592b..381877b71 100755 --- a/CUBIT_GEOCUBIT/geocubitlib/local_volume.py +++ b/CUBIT_GEOCUBIT/geocubitlib/local_volume.py @@ -22,6 +22,8 @@ # 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. # # # ############################################################################# +from __future__ import print_function + try: import start as start cubit = start.start_cubit() @@ -29,7 +31,7 @@ try: import cubit except: - print 'error importing cubit, check if cubit is installed' + print('error importing cubit, check if cubit is installed') pass numpy = start.start_numpy() @@ -39,7 +41,7 @@ def check_orientation(grdfileNAME): try: grdfile = open(grdfileNAME, 'r') - print 'reading ', grdfileNAME + print('reading ', grdfileNAME) except: txt = 'check_orintation ->error reading: ' + str(grdfile) raise Exception(txt) @@ -113,7 +115,7 @@ def process_surfacefiles(iproc, nx, ny, nstep, grdfile, unit, lat_orientation): try: grdfile = open(grdfile, 'r') - # print 'reading ',grdfile + # print('reading ',grdfile) except: txt = 'error reading: ' + str(grdfile) raise Exception(txt) @@ -123,7 +125,7 @@ def process_surfacefiles(iproc, nx, ny, nstep, grdfile, unit, lat_orientation): else: rangey = range(ny - 1, -1, -1) lat_orientation = 'NORTH2SOUTH' - print lat_orientation + print(lat_orientation) for iy in rangey: for ix in range(0, nx): txt = grdfile.readline() @@ -139,18 +141,18 @@ def process_surfacefiles(iproc, nx, ny, nstep, grdfile, unit, lat_orientation): coordy[jx, jy] = y_current elev[jx, jy] = z except: - print 'error reading point ', iy * nx + ix, txt, \ - grdfile.name, ' proc ' + print('error reading point ', iy * nx + ix, txt, \ + grdfile.name, ' proc ') raise NameError('error reading point') if (nx) * (ny) != icoord: if iproc == 0: - print 'error in the surface file ' + grdfile.name + print('error in the surface file ' + grdfile.name) if iproc == 0: - print 'x points ' + str(nx) + ' y points ' + str(ny) + \ - ' tot points ' + str((nx) * (ny)) + print('x points ' + str(nx) + ' y points ' + str(ny) + \ + ' tot points ' + str((nx) * (ny))) if iproc == 0: - print 'points read in ' + grdfile.name + ': ' + str(icoord) + print('points read in ' + grdfile.name + ': ' + str(icoord)) raise NameError grdfile.close() @@ -238,7 +240,7 @@ def read_grid(filename=None): elev[:, :, inz] = cfg.depth_top else: grdfile = cfg.filename[inz - bottomsurface] - print 'reading ', cfg.filename[inz - bottomsurface] + print('reading ', cfg.filename[inz - bottomsurface]) if cfg.irregulargridded_surf: coordx, coordy, elev_1 = process_irregular_surfacefiles( iproc, nx, ny, cfg.xmin, cfg.xmax, cfg.ymin, cfg.ymax, @@ -250,11 +252,11 @@ def read_grid(filename=None): elev[:, :, inz] = elev_1[:, :] if cfg.subduction: - print 'subduction' + print('subduction') top = elev[:, :, inz] slab = elev[:, :, inz - 1] subcrit = numpy.abs(top - slab) < cfg.subduction_thres top[subcrit] = slab[subcrit] + cfg.subduction_thres - print len(top[subcrit]) + print(len(top[subcrit])) elev[:, :, inz] = top return coordx, coordy, elev, nx, ny diff --git a/CUBIT_GEOCUBIT/geocubitlib/menu.py b/CUBIT_GEOCUBIT/geocubitlib/menu.py index 5c3b74872..46631ca29 100755 --- a/CUBIT_GEOCUBIT/geocubitlib/menu.py +++ b/CUBIT_GEOCUBIT/geocubitlib/menu.py @@ -22,9 +22,11 @@ # 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. # # # ############################################################################# +from __future__ import print_function import getopt import sys + from utilities import get_cubit_version @@ -90,10 +92,10 @@ def usage(): --SEMoutput = [YourOutputDir]) """ - print (txt) + print(txt) -# print 'reading options....' +# print('reading options....') try: if hasattr(sys, 'argv'): opts, args = getopt.getopt(sys.argv[1:], "sjmohbp1", @@ -115,8 +117,8 @@ def usage(): "build_volume", "merge1", "merge2", "merge", "collect", "meshfiles="]) -except Exception, e: - print e +except Exception as e: + print(e) sys.exit() output = 'totalmesh_merged' @@ -176,7 +178,7 @@ def usage(): for o, value in opts: if o in ('--starting_tolerance'): starting_tolerance = float(value) - # print o, value + # print(o, value) if o in ('--save_cubfile'): save_cubfile = True if o in ('--step_tolerance'): @@ -249,22 +251,22 @@ def usage(): import start as start mpiflag, iproc, numproc, mpi = start.start_mpi() if mpiflag: - print '--------, MPI ON, parallel mesher ready' + print('--------, MPI ON, parallel mesher ready') else: - print '--------, MPI OFF, serial mesher ready' + print('--------, MPI OFF, serial mesher ready') numpy = start.start_numpy() - print '--------, Numpy ON' + print('--------, Numpy ON') cubit = start.start_cubit() - print '--------, CUBIT ON' + print('--------, CUBIT ON') sys.exit() if o in ("--cfg"): cfg_name = value try: if open(cfg_name): pass - except IOError, e: - print 'error opening ', cfg_name - print e + except IOError as e: + print('error opening ', cfg_name) + print(e) import sys sys.exit() if o == ('--surface'): @@ -288,9 +290,9 @@ def usage(): import glob nf = glob.glob(cubfiles) if len(nf) > 0: - print 'cubfiles ', nf + print('cubfiles ', nf) else: - print 'files not found: ', cubfiles + print('files not found: ', cubfiles) import sys sys.exit() if o in ("--exofiles"): @@ -342,7 +344,7 @@ def usage(): pass if o in ("--addsea"): add_sea = True - print cpuxmax, cpuymax + print(cpuxmax, cpuymax) if cpuymax: pass elif cpuy > 1: @@ -361,15 +363,15 @@ def usage(): cpuxmax = cpux else: cpuxmax = 1 - print cpuxmax, cpuymax + print(cpuxmax, cpuymax) if cpml: if not cpml_size: - print 'specify the size of the cpml boundaries' + print('specify the size of the cpml boundaries') import sys sys.exit() elif cpml_size <= 0: - print 'cpml size negative, please check the parameters' + print('cpml size negative, please check the parameters') import sys sys.exit() @@ -384,7 +386,7 @@ def usage(): txt = str(k) + ' -----> ' + \ str(d[k]) txt = txt.replace("'", "").replace('"', '') - print txt + print(txt) else: try: import start as start @@ -403,7 +405,7 @@ def usage(): except: pass elif opts == []: - print __name__ + print(__name__) usage() import sys sys.exit() diff --git a/CUBIT_GEOCUBIT/geocubitlib/mesh_volume.py b/CUBIT_GEOCUBIT/geocubitlib/mesh_volume.py index a22492ce9..8aaf2db69 100755 --- a/CUBIT_GEOCUBIT/geocubitlib/mesh_volume.py +++ b/CUBIT_GEOCUBIT/geocubitlib/mesh_volume.py @@ -22,6 +22,7 @@ # 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. # # # ############################################################################# +from __future__ import print_function try: import start as start @@ -30,7 +31,7 @@ try: import cubit except: - print "error importing cubit, check if cubit is installed" + print("error importing cubit, check if cubit is installed") pass @@ -43,8 +44,8 @@ def mesh(filename=None): if cfg.map_meshing_type == 'regularmap': mesh_layercake_regularmap(filename=filename) else: - print 'error: map_meshing_type ', cfg.map_meshing_type, \ - ' not implemented' + print('error: map_meshing_type ', cfg.map_meshing_type, \ + ' not implemented') class cubitvolume: @@ -88,8 +89,8 @@ def mesh_layercake_regularmap(filename=None): command = 'composite create curve all' cubit.cmd(command) - print '###"No valid composites can be created from the specified curves." \ - is NOT a critical ERROR.' + print('###"No valid composites can be created from the specified curves." \ + is NOT a critical ERROR.') # command = "compress all" cubit.cmd(command) @@ -121,7 +122,7 @@ def mesh_layercake_regularmap(filename=None): # interval assignement surf_or, surf_vertical, list_curve_or, list_curve_vertical, \ bottom, top = get_v_h_list(list_vol, chktop=cfg.chktop) - print 'vertical surfaces: ', surf_vertical + print('vertical surfaces: ', surf_vertical) for k in surf_vertical: command = "surface " + str(k) + " scheme submap" @@ -146,7 +147,7 @@ def mesh_layercake_regularmap(filename=None): # cubit_error_stop(iproc,command,ner) if max(ucurve_interval.values()) != min(ucurve_interval.values()): schemepave = True - print 'mesh scheme is set to pave' + print('mesh scheme is set to pave') for sk in surf_or: command = "surface " + str(sk) + " scheme pave" cubit.cmd(command) @@ -164,7 +165,7 @@ def mesh_layercake_regularmap(filename=None): # cubit_error_stop(iproc,command,ner) if max(vcurve_interval.values()) != min(vcurve_interval.values()): - print 'mesh scheme is set to pave' + print('mesh scheme is set to pave') schemepave = True for sk in surf_or: command = "surface " + str(sk) + " scheme pave" @@ -247,14 +248,14 @@ def mesh_layercake_regularmap(filename=None): # # smoothing - print iproc, 'untangling...' + print(iproc, 'untangling...') cmd = "volume all smooth scheme untangle beta 0.02 cpu 10" cubit.cmd(cmd) cmd = "smooth volume all" cubit.cmd(cmd) if cfg.smoothing: - print 'smoothing .... ' + str(cfg.smoothing) + print('smoothing .... ' + str(cfg.smoothing)) cubitcommand = 'surf all smooth scheme laplacian ' cubit.cmd(cubitcommand) cubitcommand = 'smooth surf all' @@ -283,7 +284,7 @@ def mesh_layercake_regularmap(filename=None): refinement(nvol, vol, filename=filename) # # top layer vertical coarsening - print 'coarsening top layer... ', cfg.coarsening_top_layer + print('coarsening top layer... ', cfg.coarsening_top_layer) if cfg.coarsening_top_layer: cubitcommand = 'del mesh vol ' + str(vol[-1].ID) + ' propagate' cubit.cmd(cubitcommand) @@ -316,19 +317,19 @@ def mesh_layercake_regularmap(filename=None): # import boundary_definition entities = ['face'] - print iproc, 'hex block definition...' + print(iproc, 'hex block definition...') boundary_definition.define_bc(entities, parallel=True, cpux=cfg.cpux, cpuy=cfg.cpuy, cpuxmin=0, cpuymin=0, optionsea=False) # save mesh - print iproc, 'untangling...' + print(iproc, 'untangling...') cmd = "volume all smooth scheme untangle beta 0.02 cpu 10" cubit.cmd(cmd) cmd = "smooth volume all" cubit.cmd(cmd) - print iproc, 'saving...' + print(iproc, 'saving...') savemesh(mpiflag, iproc=iproc, filename=filename) # @@ -364,8 +365,8 @@ def refinement(nvol, vol, filename=None): else: for ir in cfg.tripl: if ir == 1: - print ('interface = 1 means that the refinement' - 'interface is at the bottom of the volume') + print('interface = 1 means that the refinement' + 'interface is at the bottom of the volume') txt = ' all ' idepth = 1 cubitcommand = 'refine hex in vol ' + txt @@ -512,7 +513,7 @@ def refine_inside_curve(curves, ntimes=1, depth=1, block=1, surface=False): group_id_1 = cubit.get_id_from_name("negativejac") n1 = cubit.get_group_nodes(group_id_1) if len(n1) != 0: - print 'error, negative jacobian after the refining' + print('error, negative jacobian after the refining') def get_uv_curve(list_curve_or): diff --git a/CUBIT_GEOCUBIT/geocubitlib/quality_log.py b/CUBIT_GEOCUBIT/geocubitlib/quality_log.py index 462612d69..2e6daa557 100755 --- a/CUBIT_GEOCUBIT/geocubitlib/quality_log.py +++ b/CUBIT_GEOCUBIT/geocubitlib/quality_log.py @@ -22,6 +22,8 @@ # 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. # # # ############################################################################# +from __future__ import print_function + try: import start as start cubit = start.start_cubit() @@ -29,7 +31,7 @@ try: import cubit except: - print 'error importing cubit, check if cubit is installed' + print('error importing cubit, check if cubit is installed') pass @@ -199,6 +201,6 @@ def quality_log(tqfile=None): if toclose: totstat_file.close() - print 'max specfem3d skewness: ', mesh.max_skewness - print 'min edge length: ', mesh.min_edge_length + print('max specfem3d skewness: ', mesh.max_skewness) + print('min edge length: ', mesh.min_edge_length) return mesh.max_skewness, mesh.min_edge_length diff --git a/CUBIT_GEOCUBIT/geocubitlib/read_parameter_cfg.py b/CUBIT_GEOCUBIT/geocubitlib/read_parameter_cfg.py index 164fca2c9..3a4574a96 100755 --- a/CUBIT_GEOCUBIT/geocubitlib/read_parameter_cfg.py +++ b/CUBIT_GEOCUBIT/geocubitlib/read_parameter_cfg.py @@ -22,6 +22,8 @@ # 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. # # # ############################################################################# +from __future__ import print_function + import os @@ -45,7 +47,7 @@ def readcfg(filename=None, importmenu=False, mpiflag=False): menusurface = False single = False else: - print 'error: no configuration file' + print('error: no configuration file') import sys sys.exit() # @@ -112,16 +114,16 @@ def __str__(self): for name, value in self.items(): names.append(name) values.append(value) - print names, values + print(names, values) a = zip(names, values) a.sort() arc = '' for o in a: if o[0][0] != arc: print - print o[0], ' -> ', o[1] + print(o[0], ' -> ', o[1]) arc = o[0][0] - print __name__ + print(__name__) return '____' # @@ -209,7 +211,7 @@ def __str__(self): except: pass - # print dcfg + # print(dcfg) if dcfg['nsurf']: surface_name = [] @@ -286,7 +288,7 @@ def __str__(self): pass if dcfg['irregulargridded_surf']: - print 'test' + print('test') dcfg['xmin'] = dcfg['longitude_min'] dcfg['ymin'] = dcfg['latitude_min'] dcfg['xmax'] = dcfg['longitude_max'] diff --git a/CUBIT_GEOCUBIT/geocubitlib/sea.py b/CUBIT_GEOCUBIT/geocubitlib/sea.py index 35fa5872c..65480c9a9 100755 --- a/CUBIT_GEOCUBIT/geocubitlib/sea.py +++ b/CUBIT_GEOCUBIT/geocubitlib/sea.py @@ -22,6 +22,7 @@ # 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. # # # ############################################################################# +from __future__ import print_function def adjust_sea_layers(zvertex,sealevel,bathymetry,cfg): if cfg.seaup: diff --git a/CUBIT_GEOCUBIT/geocubitlib/start.py b/CUBIT_GEOCUBIT/geocubitlib/start.py index 8d7167929..61a4e92a9 100755 --- a/CUBIT_GEOCUBIT/geocubitlib/start.py +++ b/CUBIT_GEOCUBIT/geocubitlib/start.py @@ -22,11 +22,9 @@ # 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. # # # ############################################################################# -# -# -# -# method to call the library +from __future__ import print_function +# method to call the library def start_mpi(): """ @@ -103,7 +101,7 @@ def start_cubit(init=False): import utilities cubit.init([""]) except: - print 'error importing cubit' + print('error importing cubit') sys.exit() try: if init: @@ -146,14 +144,14 @@ def start_cubit(init=False): cubit.cmd("set echo " + cfg.echo_info) cubit.cmd("set info " + cfg.cubit_info) if version_cubit > 13 and version_cubit < 15: - print 'VERSION CUBIT ', version_cubit - print 'VERSIONs of CUBIT > 13 have bugs with merge node' - print 'the merge option is not operative with this version' - print 'please download CUBIT 15+' + print('VERSION CUBIT ', version_cubit) + print('VERSIONs of CUBIT > 13 have bugs with merge node') + print('the merge option is not operative with this version') + print('please download CUBIT 15+') else: - print 'VERSION CUBIT ', version_cubit + print('VERSION CUBIT ', version_cubit) except: - print 'error start cubit' + print('error start cubit') sys.exit() return cubit @@ -193,7 +191,7 @@ def start_numpy(): try: import numpy except: - print 'error importing numpy' - print 'please check if numpy is correctly installed' - sys.exit() + print('error importing numpy') + print('please check if numpy is correctly installed') + sys.exit(1) return numpy diff --git a/CUBIT_GEOCUBIT/geocubitlib/surfaces.py b/CUBIT_GEOCUBIT/geocubitlib/surfaces.py index 872c9c459..381b70b35 100755 --- a/CUBIT_GEOCUBIT/geocubitlib/surfaces.py +++ b/CUBIT_GEOCUBIT/geocubitlib/surfaces.py @@ -22,6 +22,8 @@ # 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. # # # ############################################################################# +from __future__ import print_function + try: import start as start cubit = start.start_cubit() @@ -29,7 +31,7 @@ try: import cubit except: - print 'error importing cubit, check if cubit is installed' + print('error importing cubit, check if cubit is installed') pass @@ -76,7 +78,7 @@ def create_line_u(ind, n, step, data, unit): x, y = geo2utm(lon, lat, unit) txt = ' Position ' + str(x) + ' ' + str(y) + ' ' + str(z) command = command + txt - # print command + # print(command) cubit.silent_cmd(command) last_curve = cubit.get_last_id("curve") if last_curve != last_curve_store: @@ -94,7 +96,7 @@ def create_line_v(ind, n, n2, step, data, unit): x, y = geo2utm(lon, lat, unit) txt = ' Position ' + str(x) + ' ' + str(y) + ' ' + str(z) command = command + txt - # print command + # print(command) cubit.silent_cmd(command) last_curve = cubit.get_last_id("curve") if last_curve != last_curve_store: @@ -263,7 +265,7 @@ def create_line(position): n = index umax = max(curveskin) umin = min(curveskin) - print 'create surface skin curve ' + str(umin) + ' to ' + str(umax) + print('create surface skin curve ' + str(umin) + ' to ' + str(umax)) cubitcommand = 'create surface skin curve ' + \ str(umin) + ' to ' + str(umax) cubit.cmd(cubitcommand) diff --git a/CUBIT_GEOCUBIT/geocubitlib/utilities.py b/CUBIT_GEOCUBIT/geocubitlib/utilities.py index e4411ca9e..7eaf51f4c 100755 --- a/CUBIT_GEOCUBIT/geocubitlib/utilities.py +++ b/CUBIT_GEOCUBIT/geocubitlib/utilities.py @@ -22,6 +22,8 @@ # 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. # # # ############################################################################# +from __future__ import print_function + try: import start as start cubit = start.start_cubit() @@ -29,7 +31,7 @@ try: import cubit except: - print 'error importing cubit, check if cubit is installed' + print('error importing cubit, check if cubit is installed') pass @@ -80,7 +82,7 @@ def cubit_command_check(iproc, command, stop=True): flag = True er = cubit.get_error_count() cubit.cmd(command) - print command + print(command) ner = cubit.get_error_count() if ner > er: text = '"Proc: ' + str(iproc) + ' ERROR ' + str(command) + \ @@ -151,7 +153,7 @@ def runsave(meshfile, iproc, filename=None): saving = False if ind > len(total_saved) + 10: saving = False - print sum(total_saved), '/', len(total_saved), ' saved' + print(sum(total_saved), '/', len(total_saved), ' saved') info_total_saved = mpi.allgather(infosave) if isinstance(info_total_saved, int): @@ -211,7 +213,7 @@ def runsave(meshfile, iproc, filename=None): totstat_file.write(str(total_max_skew)) totstat_file.close() - print 'meshing process end... proc ', iproc + print('meshing process end... proc ', iproc) def importgeometry(geometryfile, iproc=0, filename=None): @@ -220,13 +222,13 @@ def importgeometry(geometryfile, iproc=0, filename=None): mpiflag, iproc, numproc, mpi = start.start_mpi() if iproc == 0: - print 'importing geometry....' + print('importing geometry....') a = ['ok from ' + str(iproc)] mpi.barrier() total_a = mpi.allgather(a) if iproc == 0: - print total_a + print(total_a) def runimport(geometryfile, iproc, filename=None): import start as start @@ -263,7 +265,7 @@ def load_curves(acis_filename): import os # # - print acis_filename + print(acis_filename) if acis_filename and os.path.exists(acis_filename): tmp_curve = cubit.get_last_id("curve") command = "import acis '" + acis_filename + "'" @@ -272,7 +274,7 @@ def load_curves(acis_filename): curves = ' '.join(str(x) for x in range(tmp_curve + 1, tmp_curve_after + 1)) elif not os.path.exists(acis_filename): - print str(acis_filename) + ' not found' + print(str(acis_filename) + ' not found') curves = None return [curves] @@ -392,7 +394,7 @@ def runsave(geometryfile, iproc, filename=None): saving = False if ind > len(total_saved) + 10: saving = False - print sum(total_saved), '/', len(total_saved), ' saved' + print(sum(total_saved), '/', len(total_saved), ' saved') info_total_saved = mpi.allgather(infosave) if isinstance(info_total_saved, int): @@ -498,10 +500,10 @@ def get_v_h_list(vol_id_list, chktop=False): # check that all the surf are Horizontal or vertical surf_all = surf_vertical + surf_or if len(surf_all) != len(lsurf): - print 'not all the surf are horizontal or vertical, check the normals' - print 'list of surfaces: ', surf_all - print 'list of vertical surface', surf_vertical - print 'list of horizontal surface', surf_or + print('not all the surf are horizontal or vertical, check the normals') + print('list of surfaces: ', surf_all) + print('list of vertical surface', surf_vertical) + print('list of horizontal surface', surf_or) bottom = [bottom] return surf_or, surf_vertical, list_curve_or, \ diff --git a/CUBIT_GEOCUBIT/geocubitlib/volumes.py b/CUBIT_GEOCUBIT/geocubitlib/volumes.py index 6bd875297..9bb50402b 100755 --- a/CUBIT_GEOCUBIT/geocubitlib/volumes.py +++ b/CUBIT_GEOCUBIT/geocubitlib/volumes.py @@ -22,6 +22,8 @@ # 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. # # # ############################################################################# +from __future__ import print_function + try: import start as start cubit = start.start_cubit() @@ -29,16 +31,16 @@ try: import cubit except: - print 'error importing cubit, check if cubit is installed' + print('error importing cubit, check if cubit is installed') pass def volumes(filename=None): """create the volumes""" import start as start - print'volume' + print('volume') cfg = start.start_cfg(filename=filename) - # print cfg + # print(cfg) # sandwich = 'verticalsandwich_volume_ascii_regulargrid_mpiregularmap' if cfg.volume_type == 'layercake_volume_ascii_regulargrid_regularmap': @@ -79,9 +81,9 @@ def surfaces(filename=None,): """create the volumes""" import start as start sandwich = 'verticalsandwich_volume_ascii_regulargrid_mpiregularmap' - print'volume' + print('volume') cfg = start.start_cfg(filename=filename) - print cfg + print(cfg) if cfg.volume_type == 'layercake_volume_ascii_regulargrid_regularmap': layercake_volume_ascii_regulargrid_mpiregularmap( filename=filename, onlysurface=True) @@ -134,7 +136,7 @@ def layercake_volume_ascii_regulargrid_mpiregularmap(filename=None, if iproc == 0 or not mpiflag: coordx_0, coordy_0, elev_0, nx_0, ny_0 = lvolume.read_grid( filename) - print 'end reading grd files ' + str(nx_0 * ny_0) + ' points' + print('end reading grd files ' + str(nx_0 * ny_0) + ' points') else: pass if iproc == 0 or not mpiflag: @@ -159,7 +161,7 @@ def layercake_volume_ascii_regulargrid_mpiregularmap(filename=None, ny = mpi.bcast() else: coordx, coordy, elev, nx, ny = lvolume.read_grid(filename) - print str(iproc) + ' end of receving grd files ' + print(str(iproc) + ' end of receving grd files ') nx_segment = int(nx / cfg.nproc_xi) + 1 ny_segment = int(ny / cfg.nproc_eta) + 1 @@ -167,13 +169,13 @@ def layercake_volume_ascii_regulargrid_mpiregularmap(filename=None, if cfg.depth_bottom != cfg.zdepth[0]: if iproc == 0: - print 'the bottom of the block is at different \ - depth than depth[0] in the configuration file' + print('the bottom of the block is at different \ + depth than depth[0] in the configuration file') nx = cfg.nproc_xi + 1 ny = cfg.nproc_eta + 1 nx_segment = 2 ny_segment = 2 - # if iproc == 0: print nx,ny,cfg.cpux,cfg.cpuy + # if iproc == 0: print(nx,ny,cfg.cpux,cfg.cpuy) # xp = (cfg.xmax - cfg.xmin) / float((nx - 1)) # yp = (cfg.ymax - cfg.ymin) / float((ny - 1)) # @@ -209,10 +211,10 @@ def layercake_volume_ascii_regulargrid_mpiregularmap(filename=None, coordx[ix, iy] = cfg.xmin + xlength * (ix) coordy[ix, iy] = cfg.ymin + ylength * (iy) - # print coordx,coordy,nx,ny + # print(coordx,coordy,nx,ny) # - print 'end of building grid ' + str(iproc) - print 'number of point: ', len(coordx) * len(coordy) + print('end of building grid ' + str(iproc)) + print('number of point: ', len(coordx) * len(coordy)) # # !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! # for each processor @@ -222,11 +224,11 @@ def layercake_volume_ascii_regulargrid_mpiregularmap(filename=None, nxmax_cpu = min(nx - 1, (nx_segment - 1) * (icpux + 1)) nymax_cpu = min(ny - 1, (ny_segment - 1) * (icpuy + 1)) # if iproc == 0: - # print nx_segment,ny_segment,nx,ny - # print icpux,icpuy,nxmin_cpu,nxmax_cpu - # print icpux,icpuy,nymin_cpu,nymax_cpu - # print coordx[0,0],coordx[nx-1,ny-1] - # print coordy[0,0],coordy[nx-1,ny-1] + # print(nx_segment,ny_segment,nx,ny) + # print(icpux,icpuy,nxmin_cpu,nxmax_cpu) + # print(icpux,icpuy,nymin_cpu,nymax_cpu) + # print(coordx[0,0],coordx[nx-1,ny-1]) + # print(coordy[0,0],coordy[nx-1,ny-1]) # # # icurve = 0 @@ -249,7 +251,7 @@ def layercake_volume_ascii_regulargrid_mpiregularmap(filename=None, else: sealevel = False bathymetry = False - print sealevel, bathymetry + print(sealevel, bathymetry) if cfg.bottomflat and inz == 0: # bottom layer # @@ -391,10 +393,10 @@ def layercake_volume_ascii_regulargrid_mpiregularmap(filename=None, # isurf = isurf + 1 else: - print "top_flat is not a valid option for sandwich" + print("top_flat is not a valid option for sandwich") else: - print inz, 'layer' + print(inz, 'layer') if cfg.geometry_format == 'regmesh': if verticalsandwich: zvertex = cfg.xwidth[inz] @@ -494,7 +496,7 @@ def layercake_volume_ascii_regulargrid_mpiregularmap(filename=None, ' ' + str(y_current) + ' ' + str(zvertex)) # - print 'proc', iproc, 'vertex list created....', len(vertex) + print('proc', iproc, 'vertex list created....', len(vertex)) uline = [] vline = [] @@ -511,7 +513,7 @@ def layercake_volume_ascii_regulargrid_mpiregularmap(filename=None, iv = iv + 1 command = 'create curve spline ' + positionx cubit.cmd(command) - # print command + # print(command) uline.append(cubit.get_last_id("curve")) for ix in range(0, nxmax_cpu - nxmin_cpu + 1): positiony = '' @@ -520,7 +522,7 @@ def layercake_volume_ascii_regulargrid_mpiregularmap(filename=None, vertex[ix + iy * (nxmax_cpu - nxmin_cpu + 1)] command = 'create curve spline ' + positiony cubit.cmd(command) - # print command + # print(command) vline.append(cubit.get_last_id("curve")) # cubit.cmd("set info " + cfg.cubit_info) @@ -528,8 +530,8 @@ def layercake_volume_ascii_regulargrid_mpiregularmap(filename=None, cubit.cmd("set journal " + cfg.jou_info) # # - print 'proc', iproc, 'lines created....', \ - len(uline), '*', len(vline) + print('proc', iproc, 'lines created....', \ + len(uline), '*', len(vline)) umax = max(uline) umin = min(uline) vmax = max(vline) @@ -573,10 +575,10 @@ def layercake_volume_ascii_regulargrid_mpiregularmap(filename=None, create_volume(inz, inz + 1, method='loft') ner2 = cubit.get_error_count() if ner != ner2: - print 'ERROR creating volume' + print('ERROR creating volume') break else: - print 'ERROR creating volume' + print('ERROR creating volume') break if ner == ner2 and not cfg.debug_geometry: # cubitcommand= 'del surface 1 to '+ str( cfg.nz ) @@ -694,7 +696,7 @@ def translate2original(xmin, ymin): # translate xmin, ymin = translate2zero() - print 'translate ...', -xmin, -ymin + print('translate ...', -xmin, -ymin) xmin_cpu = xmin_cpu - xmin ymin_cpu = ymin_cpu - ymin xmax_cpu = xmax_cpu - xmin @@ -702,7 +704,7 @@ def translate2original(xmin, ymin): ss = cubit.parse_cubit_list('surface', 'all') box = cubit.get_total_bounding_box("surface", ss) - print 'dimension... ', box + print('dimension... ', box) # cutting the surfaces xwebcut(xmin_cpu) xwebcut(xmax_cpu) diff --git a/EXAMPLES/CPML_examples/homogeneous_halfspace_HEX8_acoustic_absorbing_CPML_5sides/block_mesh-anisotropic.py b/EXAMPLES/CPML_examples/homogeneous_halfspace_HEX8_acoustic_absorbing_CPML_5sides/block_mesh-anisotropic.py index cc0a13a6a..4da6f936f 100755 --- a/EXAMPLES/CPML_examples/homogeneous_halfspace_HEX8_acoustic_absorbing_CPML_5sides/block_mesh-anisotropic.py +++ b/EXAMPLES/CPML_examples/homogeneous_halfspace_HEX8_acoustic_absorbing_CPML_5sides/block_mesh-anisotropic.py @@ -1,4 +1,5 @@ #!/usr/bin/env python +from __future__ import print_function import cubit import boundary_definition diff --git a/EXAMPLES/CPML_examples/homogeneous_halfspace_HEX8_acoustic_absorbing_CPML_6sides/block_mesh-anisotropic.py b/EXAMPLES/CPML_examples/homogeneous_halfspace_HEX8_acoustic_absorbing_CPML_6sides/block_mesh-anisotropic.py index cc0a13a6a..4da6f936f 100755 --- a/EXAMPLES/CPML_examples/homogeneous_halfspace_HEX8_acoustic_absorbing_CPML_6sides/block_mesh-anisotropic.py +++ b/EXAMPLES/CPML_examples/homogeneous_halfspace_HEX8_acoustic_absorbing_CPML_6sides/block_mesh-anisotropic.py @@ -1,4 +1,5 @@ #!/usr/bin/env python +from __future__ import print_function import cubit import boundary_definition diff --git a/EXAMPLES/CPML_examples/homogeneous_halfspace_HEX8_acoustic_elastic_absorbing_CPML_5sides/block_mesh-anisotropic.py b/EXAMPLES/CPML_examples/homogeneous_halfspace_HEX8_acoustic_elastic_absorbing_CPML_5sides/block_mesh-anisotropic.py index cc0a13a6a..4da6f936f 100755 --- a/EXAMPLES/CPML_examples/homogeneous_halfspace_HEX8_acoustic_elastic_absorbing_CPML_5sides/block_mesh-anisotropic.py +++ b/EXAMPLES/CPML_examples/homogeneous_halfspace_HEX8_acoustic_elastic_absorbing_CPML_5sides/block_mesh-anisotropic.py @@ -1,4 +1,5 @@ #!/usr/bin/env python +from __future__ import print_function import cubit import boundary_definition diff --git a/EXAMPLES/CPML_examples/homogeneous_halfspace_HEX8_acoustic_elastic_absorbing_CPML_6sides/block_mesh-anisotropic.py b/EXAMPLES/CPML_examples/homogeneous_halfspace_HEX8_acoustic_elastic_absorbing_CPML_6sides/block_mesh-anisotropic.py index cc0a13a6a..4da6f936f 100755 --- a/EXAMPLES/CPML_examples/homogeneous_halfspace_HEX8_acoustic_elastic_absorbing_CPML_6sides/block_mesh-anisotropic.py +++ b/EXAMPLES/CPML_examples/homogeneous_halfspace_HEX8_acoustic_elastic_absorbing_CPML_6sides/block_mesh-anisotropic.py @@ -1,4 +1,5 @@ #!/usr/bin/env python +from __future__ import print_function import cubit import boundary_definition diff --git a/EXAMPLES/CPML_examples/homogeneous_halfspace_HEX8_elastic_absorbing_CPML_5sides/block_mesh-anisotropic.py b/EXAMPLES/CPML_examples/homogeneous_halfspace_HEX8_elastic_absorbing_CPML_5sides/block_mesh-anisotropic.py index cc0a13a6a..4da6f936f 100755 --- a/EXAMPLES/CPML_examples/homogeneous_halfspace_HEX8_elastic_absorbing_CPML_5sides/block_mesh-anisotropic.py +++ b/EXAMPLES/CPML_examples/homogeneous_halfspace_HEX8_elastic_absorbing_CPML_5sides/block_mesh-anisotropic.py @@ -1,4 +1,5 @@ #!/usr/bin/env python +from __future__ import print_function import cubit import boundary_definition diff --git a/EXAMPLES/CPML_examples/homogeneous_halfspace_HEX8_elastic_absorbing_CPML_6sides/block_mesh-anisotropic.py b/EXAMPLES/CPML_examples/homogeneous_halfspace_HEX8_elastic_absorbing_CPML_6sides/block_mesh-anisotropic.py index cc0a13a6a..4da6f936f 100755 --- a/EXAMPLES/CPML_examples/homogeneous_halfspace_HEX8_elastic_absorbing_CPML_6sides/block_mesh-anisotropic.py +++ b/EXAMPLES/CPML_examples/homogeneous_halfspace_HEX8_elastic_absorbing_CPML_6sides/block_mesh-anisotropic.py @@ -1,4 +1,5 @@ #!/usr/bin/env python +from __future__ import print_function import cubit import boundary_definition diff --git a/EXAMPLES/Gmsh_simple_lddrk/Gmsh2specfem.py b/EXAMPLES/Gmsh_simple_lddrk/Gmsh2specfem.py deleted file mode 100755 index 785d0a648..000000000 --- a/EXAMPLES/Gmsh_simple_lddrk/Gmsh2specfem.py +++ /dev/null @@ -1,510 +0,0 @@ -#!/usr/bin/env python -# -# uses meshio: -# https://github.com/nschloe/meshio -# -# install with: -# pip install -U meshio -# -# this script reads in a Gmsh file (in Gmsh-ASCII format) and outputs SPECFEM-readable mesh files -# -# Boundary surfaces must be defined prior when generating the mesh as Physical Surfaces in Gmsh. -# Surfaces must have names: -# - "top" -# - "bottom" -# - "xmin" -# - "xmax" -# - "ymin" -# - "ymax" -# -from __future__ import print_function -import sys -import os -import math - -import numpy as np - -try: - import meshio -except: - print("importing module meshio failed, please make sure to install it via: pip install -U meshio") - sys.exit(1) - -#-------------------------------------------------------------------------------------------------- - -def read_mesh_file_msh(file): - """ - reads in Gmsh file - """ - print("") - print("reading file: ",file) - print("") - - # reads *.msh file (Gmsh Ascii-format) - mesh = meshio.read(file) - points, cells, point_data, cell_data, field_data = ( - mesh.points, mesh.cells, - mesh.point_data, mesh.cell_data, mesh.field_data) - - print("mesh data:") - print("number of points: ",len(points)) - print("cells : ",len(cells),"items") - for name,_ in cells.items(): - print(" ",name) - - print("point_data: ",len(point_data),"items") - for name,dic in point_data.items(): - print(" ",name) - for s,array in dic.items(): - print(" ",s,len(array)) - - print("cell_data : ",len(cell_data),"items") - for name,dic in cell_data.items(): - print(" ",name) - for s,array in dic.items(): - print(" ",s,len(array)) - - print("field_data: ",len(field_data),"items") - for name,_ in field_data.items(): - print(" ",name) - - print("") - - #debug - #print("cells:") - #print(cells) - #print("cell_data:") - #print(cell_data) - #print("field_data:") - #print(field_data) - #print("") - - # cleans points array (removing unused points) - points = clean_out_unused_points(points,cells) - - return points,cells,point_data,cell_data,field_data - -#-------------------------------------------------------------------------------------------------- - - -def clean_out_unused_points(points,cells): - # for some reason, there might be extra nodes stored. we will clean them out and only use the onces needed by the hexas. - # assumes that the surface quads won't require new nodes. - print("cleaning points...") - - # checks if anything to do - if not 'hexahedron' in cells.keys(): - return points - - # get array with all point entries - cell_all = cells['hexahedron'].flatten() - - # remove redundant entries - cell_points = (list(set(cell_all))) - - # point indices [0,1,2,..] - points_index = range(len(points)) - - # array with points not listed by cells (difference between the two sets) - unused_points = list(set(points_index) - set(cell_points)) - - print(" number of points = ",len(points)) - print(" number of points needed by hexahedra = ",len(cell_points)) - print(" unused points: ",unused_points) - - # checks if anything to do - if len(unused_points) > 0: - # remove unused points from lists - points = np.delete(points, unused_points, axis=0) - - print(" new number of points = ",len(points)) - - # not used: updates point data - #for key in point_data: - # point_data[key] = numpy.delete(point_data[key],unused_points,axis=0) - - return points - - - -#-------------------------------------------------------------------------------------------------- - - -def get_surface_dimension(surface,points): - # determines min/max dimensions of all points on a surface - surf_points_x = [] - surf_points_y = [] - surf_points_z = [] - for elem in surface: - #print(elem) - for i in range(1,5): - index = elem[i] - if index < 0 or index >= len(points): - print("invalid point index ",index,"in surface") - sys.exit(1) - #print(elem,index,points[index]) - surf_points_x.append(points[index][0]) - surf_points_y.append(points[index][1]) - surf_points_z.append(points[index][2]) - - xmin = min(surf_points_x) - xmax = max(surf_points_x) - ymin = min(surf_points_y) - ymax = max(surf_points_y) - zmin = min(surf_points_z) - zmax = max(surf_points_z) - - #zmax = max(p[2] for p in surf_points) - - #debug - #print(surface) - #print(surf_points_x) - #print(surf_points_y) - #print(surf_points_z) - - return xmin,xmax,ymin,ymax,zmin,zmax - -#-------------------------------------------------------------------------------------------------- - -def detect_boundary_surfaces(points,cells,cell_data,field_data): - """ - determines surfaces at xmin,xmax,ymin,ymax,top and bottom - """ - print("determining boundary surfaces:") - - surf_top = [] - surf_bottom = [] - surf_xmin = [] - surf_xmax = [] - surf_ymin = [] - surf_ymax = [] - - # checks - if not 'quad' in cells.keys(): - print("invalid mesh: surfaces need 'quad' cells") - sys.exit(1) - if not 'hexahedron' in cells.keys(): - print("invalid mesh: surfaces need 'hexahedron' cells") - sys.exit(1) - - quads = cells['quad'] - hexas = cells['hexahedron'] - - # checks if valid entries - if len(quads) == 0: - print("Error: need entries in 'quad' cells") - sys.exit(1) - if len(hexas) == 0: - print("Error: need entries in 'hexahedron' cells") - sys.exit(1) - - # gets surface index - index = [0,0,0,0,0,0] # format top,bottom,xmin,xmax,ymin,ymax - for key in field_data.keys(): - # top surface - if key == 'top': - # for example item: u'top': array([2, 2]) # first index is surface id - index[0] = field_data[key][0] - # bottom - if key == 'bottom': - index[1] = field_data[key][0] - # xmin - if key == 'xmin': - index[2] = field_data[key][0] - # xmax - if key == 'xmax': - index[3] = field_data[key][0] - # ymin - if key == 'ymin': - index[4] = field_data[key][0] - # ymax - if key == 'ymax': - index[5] = field_data[key][0] - print(" surface index: ",index) - print("") - - # gets cell data - if not 'quad' in cell_data.keys(): - print("invalid mesh: needs 'quad' entry in cell_data") - sys.exit(1) - quad_dict = cell_data['quad'] - - # meshio creates 'gmsh:physical' and 'gmsh:geometrical' items, using physical volume associations - if not 'gmsh:physical' in quad_dict.keys(): - print("invalid mesh: needs 'gmsh:physical' entry in cell_data for 'quads'") - sys.exit(1) - quad_mat = quad_dict['gmsh:physical'] - print(" number of quads with material = ",len(quad_mat)) - print("") - - if len(quad_mat) == 0: - print("Error: need quads with material for boundary surfaces, please define physical surfaces using Gmsh") - sys.exit(1) - if len(quad_mat) != len(quads): - print("invalid total number of quads ",len(quads),"and number of quads with cell data ",len(quad_mat)) - sys.exit(1) - - index_count = [0,0,0,0,0,0] - for i,q in enumerate(quad_mat): - ind = quad_mat[i] - if ind == index[0]: index_count[0] += 1 - if ind == index[1]: index_count[1] += 1 - if ind == index[2]: index_count[2] += 1 - if ind == index[3]: index_count[3] += 1 - if ind == index[4]: index_count[4] += 1 - if ind == index[5]: index_count[5] += 1 - print(" surface index counts: ",index_count) - print("") - - # the following assumes linear elements, i.e., quads defined by 4 corners only - if len(quads[0]) != 4: - print("Invalid element type for quads: need linear elements, higher order elements not supported yet") - - # loops through quads - for i,q in enumerate(quads): - #print(i,q) # i starts at 0,1,2,..; quads[id1,id2,id3,id4] indices also start at 0 - - # for specfem, we will need to find corresponding hexahedral element which contains this quad - # DK DK significantly faster way suggested by Deyu Ming - #el = [ [j,h] for j,h in enumerate(hexas) if q[0] in h and q[1] in h and q[2] in h and q[3] in h ] - el = np.where(np.count_nonzero(np.isin(hexas,q),axis=1)==4) - - # boundary surfaces should only associate with a single hexahedral element - if len(el) != 1: - print("Error: could not find hexahedral for surface quad ",i,q,"el = ",el) - sys.exit(1) - hex_index = el[0][0] - - # adds surface id to quad corners - elem = [hex_index,q[0],q[1],q[2],q[3]] - - # gets associated surface type index - ind = quad_mat[i] - - # appends to surface - if ind == index[0]: surf_top.append(elem) - if ind == index[1]: surf_bottom.append(elem) - if ind == index[2]: surf_xmin.append(elem) - if ind == index[3]: surf_xmax.append(elem) - if ind == index[4]: surf_ymin.append(elem) - if ind == index[5]: surf_ymax.append(elem) - - print(" top surface: number of quads = ",len(surf_top)) - if len(surf_top) > 0: - xmin,xmax,ymin,ymax,zmin,zmax = get_surface_dimension(surf_top,points) - print(" dimension:") - print(" xmin/xmax = ",xmin,"/",xmax) - print(" ymin/ymax = ",ymin,"/",ymax) - print(" zmin/zmax = ",zmin,"/",zmax) - - print(" bottom surface: number of quads = ",len(surf_bottom)) - if len(surf_bottom) > 0: - xmin,xmax,ymin,ymax,zmin,zmax = get_surface_dimension(surf_bottom,points) - print(" dimension:") - print(" xmin/xmax = ",xmin,"/",xmax) - print(" ymin/ymax = ",ymin,"/",ymax) - print(" zmin/zmax = ",zmin,"/",zmax) - - print(" xmin surface: number of quads = ",len(surf_xmin)) - if len(surf_xmin) > 0: - xmin,xmax,ymin,ymax,zmin,zmax = get_surface_dimension(surf_xmin,points) - print(" dimension:") - print(" xmin/xmax = ",xmin,"/",xmax) - print(" ymin/ymax = ",ymin,"/",ymax) - print(" zmin/zmax = ",zmin,"/",zmax) - - print(" xmax surface: number of quads = ",len(surf_xmax)) - if len(surf_xmax) > 0: - xmin,xmax,ymin,ymax,zmin,zmax = get_surface_dimension(surf_xmax,points) - print(" dimension:") - print(" xmin/xmax = ",xmin,"/",xmax) - print(" ymin/ymax = ",ymin,"/",ymax) - print(" zmin/zmax = ",zmin,"/",zmax) - - print(" ymin surface: number of quads = ",len(surf_ymin)) - if len(surf_ymin) > 0: - xmin,xmax,ymin,ymax,zmin,zmax = get_surface_dimension(surf_ymin,points) - print(" dimension:") - print(" xmin/xmax = ",xmin,"/",xmax) - print(" ymin/ymax = ",ymin,"/",ymax) - print(" zmin/zmax = ",zmin,"/",zmax) - - print(" ymax surface: number of quads = ",len(surf_ymax)) - if len(surf_ymax) > 0: - xmin,xmax,ymin,ymax,zmin,zmax = get_surface_dimension(surf_ymax,points) - print(" dimension:") - print(" xmin/xmax = ",xmin,"/",xmax) - print(" ymin/ymax = ",ymin,"/",ymax) - print(" zmin/zmax = ",zmin,"/",zmax) - - print("") - - return surf_top,surf_bottom,surf_xmin,surf_xmax,surf_ymin,surf_ymax - -#-------------------------------------------------------------------------------------------------- - -def export_mesh(points,cells,point_data,cell_data,field_data): - """ - creates files in MESH/ directory - """ - print("") - print("exporting mesh files:") - - # outputs mesh files - os.system("mkdir -p MESH/") - - # node coordinates - print(" number of nodes = ",len(points)) - filename = "MESH/nodes_coords_file" - with open(filename,"w") as f: - f.write("%d\n" % len(points)) - for i,p in enumerate(points): - # format: #id_node #x_coordinate #y_coordinate #z_coordinate - f.write("%d %f %f %f\n" % (i+1,p[0],p[1],p[2]) ) - print(" exported to: ",filename) - print("") - - # cells - if not 'hexahedron' in cells.keys(): - print("invalid mesh: needs 'hexahedron' cells") - sys.exit(1) - hexas = cells['hexahedron'] - print(" number of hexahedra = ",len(hexas)) - filename = "MESH/mesh_file" - with open(filename,"w") as f: - f.write("%d\n" % len(hexas)) - for i,h in enumerate(hexas): - # format: # element_id #id_node1 ... #id_node8 - # note: specfem assumes indices start at 1 (fortran style, not at 0 like in python or C) - f.write("%d %d %d %d %d %d %d %d %d\n" % (i+1,h[0]+1,h[1]+1,h[2]+1,h[3]+1,h[4]+1,h[5]+1,h[6]+1,h[7]+1) ) - print(" exported to: ",filename) - print("") - - # material (or volume) associations - if not 'hexahedron' in cell_data.keys(): - print("invalid mesh: needs 'hexahedron' entry in cell_data") - sys.exit(1) - hex_dict = cell_data['hexahedron'] - - # meshio creates 'gmsh:physical' and 'gmsh:geometrical' items, using physical volume associations - if not 'gmsh:physical' in hex_dict.keys(): - print("invalid mesh: needs 'gmsh:physical' entry in cell_data for 'hexahedron'") - sys.exit(1) - hex_mat = hex_dict['gmsh:physical'] - print(" number of hexahedra with material = ",len(hex_mat)) - # checks if number of physical hexahedras match - if len(hex_mat) != len(hexas): - print("invalid numbers of hexahedras ",len(hexas),"and material associations ",len(hex_mat)) - print("please define physical volumes using Gmsh") - sys.exit(1) - filename = "MESH/materials_file" - with open(filename,"w") as f: - for i,mat in enumerate(hex_mat): - # format: #id_element #flag - f.write("%d %d\n" % (i+1,mat) ) - print(" exported to: ",filename) - print("") - - # creates a default material - # format: #(1)material_domain_id #(2)material_id #(3)rho #(4)vp #(5)vs #(6)Q_kappa #(7)Q_mu #(8)anisotropy_flag - - # default values - domain_id = 2 # 1==acoustic/2==elastic - rho = 2300.0 - vp = 2800.0 - vs = 1500.0 - q_kappa = 9999.0 - q_mu = 9999.0 - ani_flag = 0 - - # checks mesh definitions - if len(hex_mat) == 0: - print("Error: need hexahedra with material, please define physical volumes using Gmsh") - sys.exit(1) - - material_ids = [hex_mat[0]] - for mat in hex_mat: - if not mat in material_ids: material_ids.append(mat) - print(" material ids: ",material_ids) - - filename = "MESH/nummaterial_velocity_file" - if not os.path.isfile(filename): - print(" creating a default nummaterial velocity file...") - with open(filename,"w") as f: - for mat_id in material_ids: - # format: #(1)material_domain_id #(2)material_id #(3)rho #(4)vp #(5)vs #(6)Q_kappa #(7)Q_mu #(8)anisotropy_flag - f.write("%d %d %f %f %f %f %f %d\n" % (domain_id,mat_id,rho,vp,vs,q_kappa,q_mu,ani_flag)) - # adds info - f.write("\n") - f.write("! note: format of nummaterial_velocity_file must be\n") - f.write("! #(1)material_domain_id #(2)material_id #(3)rho #(4)vp #(5)vs #(6)Q_kappa #(7)Q_mu #(8)anisotropy_flag\n") - print(" exported to: ",filename) - print("") - - - # surface elements - print("surfaces:") - if not 'quad' in cells.keys(): - print("invalid mesh: surfaces need 'quad' cells") - sys.exit(1) - quads = cells['quad'] - print(" number of quads = ",len(quads)) - print("") - - surf_top,surf_bottom,surf_xmin,surf_xmax,surf_ymin,surf_ymax = detect_boundary_surfaces(points,cells,cell_data,field_data) - - # surfaces - short = ["zmax", "bottom", "xmin", "xmax", "ymin", "ymax" ] - filenames = ["MESH/free_or_absorbing_surface_file_zmax", "MESH/absorbing_surface_file_bottom", - "MESH/absorbing_surface_file_xmin", "MESH/absorbing_surface_file_xmax", - "MESH/absorbing_surface_file_ymin", "MESH/absorbing_surface_file_ymax" ] - surfaces = [ surf_top, surf_bottom, surf_xmin, surf_xmax, surf_ymin, surf_ymax ] - - for name,filename,surface in zip(short,filenames,surfaces): - # info - print("surface {:>8}: {}".format(name,len(surface))) - if len(surface) > 0: - with open(filename,"w") as f: - f.write("%d\n" % len(surface)) - for elem in surface: - # format: #id_(element containing the face) #id_node1_face .. #id_node4_face - # note: specfem assumes indices start at 1 (fortran style, not at 0 like in python or C) - f.write("%d %d %d %d %d\n" % (elem[0]+1,elem[1]+1,elem[2]+1,elem[3]+1,elem[4]+1)) - print(" exported to: ",filename) - print("") - - print("") - print(" done, see files in directory: MESH/") - print("") - return - -#-------------------------------------------------------------------------------------------------- - -def export2SPECFEM3D(filename): - """ - converts Gmsh file to SPECFEM mesh files - """ - # reads in mesh file - points,cells,point_data,cell_data,field_data = read_mesh_file_msh(filename) - - # exports in specfem format - export_mesh(points,cells,point_data,cell_data,field_data) - - -def usage(): - print("usage: ./Gmsh2specfem.py file") - print(" where") - print(" file - a Gmsh file, *.msh (Gmsh Ascii-format)") - - -if __name__ == '__main__': - # gets argument - if len(sys.argv) != 2: - usage() - sys.exit(1) - else: - filename = sys.argv[1] - - # main routine - export2SPECFEM3D(filename) - diff --git a/EXAMPLES/Gmsh_simple_lddrk/README b/EXAMPLES/Gmsh_simple_lddrk/README index 85d744166..85a8f92be 100644 --- a/EXAMPLES/Gmsh_simple_lddrk/README +++ b/EXAMPLES/Gmsh_simple_lddrk/README @@ -10,8 +10,10 @@ requires: * pygmsh - python wrapper, see https://pypi.python.org/pypi/pygmsh * numpy - scientific computing with python, see http://www.numpy.org -tested for Gmsh version 3.0, pygmsh version 4.1.0, numpy version 1.14.0 +tested for: +- python 2.7: Gmsh version 3.0, pygmsh version 4.1.0, numpy version 1.14.0 +- python 3.7: Gmsh version 4.0, pygmsh version 6.0.2, numpy version 1.17.2 0. installation: diff --git a/EXAMPLES/Gmsh_simple_lddrk/mesh_example_with_gmsh.py b/EXAMPLES/Gmsh_simple_lddrk/mesh_example_with_gmsh.py index d0d631cb1..cd39207e7 100755 --- a/EXAMPLES/Gmsh_simple_lddrk/mesh_example_with_gmsh.py +++ b/EXAMPLES/Gmsh_simple_lddrk/mesh_example_with_gmsh.py @@ -19,6 +19,7 @@ # it uses gmsh to generate a hexahedral mesh. # from __future__ import print_function + import sys import os import math @@ -34,6 +35,7 @@ import meshio # from python file import +sys.path.append('../../utils/Cubit_or_Gmsh/') from Gmsh2specfem import export2SPECFEM3D @@ -63,6 +65,9 @@ ############################################################## +# globals +python_major_version = 0 +pygmsh_major_version = 0 #-------------------------------------------------------------------------------------------------- # @@ -277,6 +282,7 @@ def mesh_3D(): # thus to obtain a top with some elevation, we will stretch the mesh points accordingly after the mesh was generated. # it is done here in a simple way to show how one could modify the mesh. global xsize,ysize,zsize,mesh_element_size_XY,mesh_element_size_Z + global python_major_version,pygmsh_major_version # output directory for mesh files os.system("mkdir -p MESH/") @@ -353,19 +359,31 @@ def mesh_3D(): geom.add_raw_code('Coherence;') # Physical Volume - geom.add_physical_volume(vol, label="vol1") + if pygmsh_major_version < 6: + geom.add_physical_volume(vol, label="vol1") + else: + geom.add_physical(vol, label="vol1") # Physical Surfaces # note: extrude only returns lateral surfaces. we thus had to create the bottom surface explicitly. # geometries and mesh from extruded surface will then merge with bottom surface. # top and bottom - geom.add_physical_surface(surface_top, label='top') - geom.add_physical_surface(surface_bottom, label='bottom') - # lateral sides - geom.add_physical_surface(lat[0], label='ymin') - geom.add_physical_surface(lat[1], label='xmax') - geom.add_physical_surface(lat[2], label='ymax') - geom.add_physical_surface(lat[3], label='xmin') + if pygmsh_major_version < 6: + geom.add_physical_surface(surface_top, label='top') + geom.add_physical_surface(surface_bottom, label='bottom') + # lateral sides + geom.add_physical_surface(lat[0], label='ymin') + geom.add_physical_surface(lat[1], label='xmax') + geom.add_physical_surface(lat[2], label='ymax') + geom.add_physical_surface(lat[3], label='xmin') + else: + geom.add_physical(surface_top, label='top') + geom.add_physical(surface_bottom, label='bottom') + # lateral sides + geom.add_physical(lat[0], label='ymin') + geom.add_physical(lat[1], label='xmax') + geom.add_physical(lat[2], label='ymax') + geom.add_physical(lat[3], label='xmin') # explicit volume meshing .. not needed, will be done when called by: gmsh -3 .. below #geom.add_raw_code('Mesh 3;') @@ -383,7 +401,11 @@ def mesh_3D(): print("") # meshing - points, cells, point_data, cell_data, field_data = pygmsh.generate_mesh(geom,verbose=True) + if pygmsh_major_version < 6: + points, cells, point_data, cell_data, field_data = pygmsh.generate_mesh(geom,verbose=True) + else: + mesh = pygmsh.generate_mesh(geom,verbose=True) + points,cells,point_data,cell_data,field_data = mesh.points,mesh.cells,mesh.point_data,mesh.cell_data,mesh.field_data # mesh info print("") @@ -431,8 +453,16 @@ def mesh_3D(): def create_mesh(): + global python_major_version,pygmsh_major_version + # version info - print("pygmsh version: ",pygmsh.__version__) + python_major_version = sys.version_info[0] + python_minor_version = sys.version_info[1] + print("Python version: ","{}.{}".format(python_major_version,python_minor_version)) + + version = pygmsh.__version__ + pygmsh_major_version = int(version.split(".")[0]) + print("pygmsh version: ",version) version = pygmsh.get_gmsh_major_version() print("Gmsh version : ",version) diff --git a/EXAMPLES/LTS_homogeneous_halfspace_HEX8/block_mesh.py b/EXAMPLES/LTS_homogeneous_halfspace_HEX8/block_mesh.py index 6f57641ae..8a9ffb141 100755 --- a/EXAMPLES/LTS_homogeneous_halfspace_HEX8/block_mesh.py +++ b/EXAMPLES/LTS_homogeneous_halfspace_HEX8/block_mesh.py @@ -1,4 +1,5 @@ #!/usr/bin/env python +from __future__ import print_function import os import sys @@ -17,7 +18,7 @@ version = cubit.get_version() version_major = int(version.split(".")[0]) version_minor = int(version.split(".")[1]) -print "cubit version: ",version +print("cubit version: ",version) cubit.cmd('reset') cubit.cmd('brick x 134000 y 134000 z 60000') diff --git a/EXAMPLES/Mount_StHelens/mesh_mount.py b/EXAMPLES/Mount_StHelens/mesh_mount.py index 40f576a16..505a3f52c 100755 --- a/EXAMPLES/Mount_StHelens/mesh_mount.py +++ b/EXAMPLES/Mount_StHelens/mesh_mount.py @@ -11,49 +11,43 @@ # # ############################################################# -import cubit -cubit.init([""]) - - -try: - -from geocubitlib import boundary_definition, exportlib -from geocubitlib import cubit2specfem3d - -except: - -import boundary_definition -import cubit2specfem3d - +from __future__ import print_function import os import sys import os.path import time +import cubit +cubit.init([""]) - +try: + from geocubitlib import boundary_definition, exportlib + from geocubitlib import cubit2specfem3d +except: + import boundary_definition + import cubit2specfem3d # time stamp -print "#" + time.strftime("%a, %d %b %Y %H:%M:%S", time.localtime()) +print("#" + time.strftime("%a, %d %b %Y %H:%M:%S", time.localtime())) # working directory cwd = os.getcwd() -print "#current working directory: " + str(cwd) +print("#current working directory: " + str(cwd)) if cwd[len(cwd)-14:len(cwd)] != "Mount_StHelens": - print "" - print "#please run this script from example directory: SPECFEM3D/example/Mount_StHelens/" - print "" + print("") + print("#please run this script from example directory: SPECFEM3D/example/Mount_StHelens/") + print("") cubit.cmd('version') cubit.cmd('reset') -print "running meshing script..." -print "" -print "note: this script uses topography surface in ACIS format" -print " meshing will take around 15 min" -print "" +print("running meshing script...") +print("") +print("note: this script uses topography surface in ACIS format") +print(" meshing will take around 15 min") +print("") # uses developer commands cubit.cmd('set developer commands on') @@ -64,18 +58,18 @@ # 0. step: loading topography surface # ############################################################# -print "#" + time.strftime("%a, %d %b %Y %H:%M:%S", time.localtime()) -print "#loading topo surface..." +print("#" + time.strftime("%a, %d %b %Y %H:%M:%S", time.localtime())) +print("#loading topo surface...") # topography surface if os.path.exists("topo.cub"): - print "opening existing topography surface" + print("opening existing topography surface") # topography surface # previously run, just reopen the cubit file cubit.cmd('open "topo.cub"') else: # topo surface doesn't exist yet, this creates it: - print "reading in topography surface" + print("reading in topography surface") # reads in topography points and creates sheet surface execfile("./read_topo.py") # clear @@ -92,8 +86,8 @@ # 1. step: creates temporary brick volume # ############################################################# -print "#" + time.strftime("%a, %d %b %Y %H:%M:%S", time.localtime()) -print "#creating brick..." +print("#" + time.strftime("%a, %d %b %Y %H:%M:%S", time.localtime())) +print("#creating brick...") # creates temporary brick volume # new brick volume (depth will become 1/2 * 20,000 = 10,000 m) @@ -110,11 +104,11 @@ # 2. step: creates volume with topography surface # ############################################################# -print "#" + time.strftime("%a, %d %b %Y %H:%M:%S", time.localtime()) -print "#creating volume with topography..." +print("#" + time.strftime("%a, %d %b %Y %H:%M:%S", time.localtime())) +print("#creating volume with topography...") -print "#" + time.strftime("%a, %d %b %Y %H:%M:%S", time.localtime()) -print "#imprinting volume, this will take around 1 min, please be patience..." +print("#" + time.strftime("%a, %d %b %Y %H:%M:%S", time.localtime())) +print("#imprinting volume, this will take around 1 min, please be patience...") cubit.cmd('merge all') cubit.cmd('imprint all') @@ -131,9 +125,9 @@ ############################################################# # checks if new file available if not os.path.exists("topo_2.acis"): - print "" - print "error creating new volume, please check manually..." - print "" + print("") + print("error creating new volume, please check manually...") + print("") cubit.cmd('pause') # clears workspace cubit.cmd('reset') @@ -141,8 +135,8 @@ # imports surfaces and merges to single volume cubit.cmd('import acis "topo_2.acis" ascii merge_globally') -print "#" + time.strftime("%a, %d %b %Y %H:%M:%S", time.localtime()) -print "#creating new volume, this will take another 2 min..." +print("#" + time.strftime("%a, %d %b %Y %H:%M:%S", time.localtime())) +print("#creating new volume, this will take another 2 min...") cubit.cmd('create volume surface all heal') # backup @@ -153,9 +147,9 @@ # 4. step: create mesh # ############################################################# -print "#" + time.strftime("%a, %d %b %Y %H:%M:%S", time.localtime()) -print "#initial meshing..." -print "#(will take around 7 min)" +print("#" + time.strftime("%a, %d %b %Y %H:%M:%S", time.localtime())) +print("#initial meshing...") +print("#(will take around 7 min)") # optional: refining mesh at surface # @@ -182,9 +176,9 @@ else: # optional surface refinement # time stamp - print "#" + time.strftime("%a, %d %b %Y %H:%M:%S", time.localtime()) - print "#refining surface mesh..." - print "#(will take around 3 hours)" + print("#" + time.strftime("%a, %d %b %Y %H:%M:%S", time.localtime())) + print("#refining surface mesh...") + print("#(will take around 3 hours)") # starts with a crude mesh elementsize = 2000.0 cubit.cmd('volume all size '+str(elementsize)) @@ -215,8 +209,8 @@ # time stamp -print "#" + time.strftime("%a, %d %b %Y %H:%M:%S", time.localtime()) -print "#done meshing..." +print("#" + time.strftime("%a, %d %b %Y %H:%M:%S", time.localtime())) +print("#done meshing...") # backup cubit.cmd('save as "topo_4.cub" overwrite') @@ -250,6 +244,6 @@ # all files needed by SCOTCH are now in directory MESH # time stamp -print "#" + time.strftime("%a, %d %b %Y %H:%M:%S", time.localtime()) +print("#" + time.strftime("%a, %d %b %Y %H:%M:%S", time.localtime())) diff --git a/EXAMPLES/Mount_StHelens/mesh_mount_stl.py b/EXAMPLES/Mount_StHelens/mesh_mount_stl.py index dec0ba8ea..16c289155 100755 --- a/EXAMPLES/Mount_StHelens/mesh_mount_stl.py +++ b/EXAMPLES/Mount_StHelens/mesh_mount_stl.py @@ -8,8 +8,15 @@ # ( try out script mesh_mount.py when using a different CUBIT version) # ############################################################# +from __future__ import print_function + +import os +import sys +import os.path + import cubit cubit.init([""]) + try: from geocubitlib import boundary_definition from geocubitlib import cubit2specfem3d @@ -17,31 +24,27 @@ import boundary_definition import cubit2specfem3d -import os -import sys -import os.path - # time stamp -print "#" + time.strftime("%a, %d %b %Y %H:%M:%S", time.localtime()) +print("#" + time.strftime("%a, %d %b %Y %H:%M:%S", time.localtime())) # working directory cwd = os.getcwd() -print "#current working directory: " + str(cwd) +print("#current working directory: " + str(cwd)) if cwd[len(cwd)-14:len(cwd)] != "Mount_StHelens": - print "" - print "#please run this script from example directory: SPECFEM3D/example/Mount_StHelens/" - print "" + print("") + print("#please run this script from example directory: SPECFEM3D/example/Mount_StHelens/") + print("") cubit.cmd('version') cubit.cmd('reset') os.system('rm -f topo_brick.stl topo_vol.stl topo_vol2.stl') -print "running meshing script..." -print "" -print "note: this script uses topography surface in STL format" -print " meshing will take around 2 min" -print "" +print("running meshing script...") +print("") +print("note: this script uses topography surface in STL format") +print(" meshing will take around 2 min") +print("") # note: this is a workaround to use STL file formats rather than ACIS formats. # for our purpose to create a simple mesh, this STL formats are faster @@ -68,9 +71,9 @@ cubit.cmd('reset') #checks if new file available if not os.path.exists("topo_brick.stl"): - print "" - print "error creating new STL file topo_brick.stl, please check manually..." - print "" + print("") + print("error creating new STL file topo_brick.stl, please check manually...") + print("") cubit.cmd('pause') ############################################################# @@ -81,12 +84,12 @@ # topography surface if os.path.exists("topo.stl"): - print "opening existing topography surface" + print("opening existing topography surface") # previously run, just reopen the cubit file cubit.cmd('import stl "topo.stl" merge stitch') else: # topo surface doesn't exist yet, this creates it: - print "reading in topography surface" + print("reading in topography surface") # reads in topography points and creates sheet surface execfile("./read_topo.py") # clear @@ -111,9 +114,9 @@ cubit.cmd('reset') #checks if new file available if not os.path.exists("topo_vol.stl"): - print "" - print "error creating new STL file topo_vol.stl, please check manually..." - print "" + print("") + print("error creating new STL file topo_vol.stl, please check manually...") + print("") cubit.cmd('pause') ############################################################# @@ -125,9 +128,9 @@ os.system('awk \'BEGIN{print \"solid Body_1\";}{if($0 !~ /solid/) print $0;}END{print \"endsolid Body_1\";}\' topo_vol.stl > topo_vol2.stl') #checks if new file available if not os.path.exists("topo_vol2.stl"): - print "" - print "error creating new STL file topo_vol2.stl, please check manually..." - print "" + print("") + print("error creating new STL file topo_vol2.stl, please check manually...") + print("") cubit.cmd('pause') ############################################################# @@ -200,4 +203,4 @@ # all files needed by SCOTCH are now in directory MESH # time stamp -print "#" + time.strftime("%a, %d %b %Y %H:%M:%S", time.localtime()) +print("#" + time.strftime("%a, %d %b %Y %H:%M:%S", time.localtime())) diff --git a/EXAMPLES/Mount_StHelens/read_topo.py b/EXAMPLES/Mount_StHelens/read_topo.py index 0c1628cae..1709d5fa0 100755 --- a/EXAMPLES/Mount_StHelens/read_topo.py +++ b/EXAMPLES/Mount_StHelens/read_topo.py @@ -1,6 +1,5 @@ #!/usr/bin/env python - -import cubit +from __future__ import print_function import os import sys @@ -8,7 +7,9 @@ import string import math -print sys.path +import cubit + +print(sys.path) ############################################################# # USER PARAMETERS @@ -25,7 +26,7 @@ #os.system('./convert_lonlat2utm.pl ptopo.mean.xyz 10 > ptopo.mean.utm ') # reads in points -print '#reading from file: ',inputFile +print('#reading from file: ',inputFile) cubit.cmd('reset') cubit.cmd('echo off') @@ -33,40 +34,40 @@ # creates point vertices -print '#creating points...' +print('#creating points...') count = 0 for line in fileinput.input( inputFile ): count = count + 1 - #print '# '+str(count)+' point: '+line + #print('# '+str(count)+' point: '+line) lineitems = line.split() x = lineitems[0] y = lineitems[1] z = lineitems[2] xyz = str(x) + ' ' + str(y) + ' ' + str(z) - #print '#point: ',xyz + #print('#point: ',xyz) cubit.cmd('create vertex '+ xyz ) fileinput.close() -print '#done points: '+str(count) -print '' +print('#done points: '+str(count)) +print('') cubit.cmd('Display') # creates smooth spline curves for surface -print '#creating curves...' +print('#creating curves...') countcurves = 0 for i in range(1,count+1): if i > 1 : if i % nstep == 0 : countcurves = countcurves + 1 cubit.cmd('create curve spline vertex '+str(i-nstep+1) + ' to ' + str(i) + ' delete' ) -print '#done curves: '+str(countcurves) -print '' +print('#done curves: '+str(countcurves)) +print('') cubit.cmd('Display') cubit.cmd('pause') # creates surface -print '#creating skin surface...' +print('#creating skin surface...') cubit.cmd('create surface skin curve all') @@ -78,7 +79,7 @@ cubit.cmd('delete vertex all') cubit.cmd('delete curve all') -print '#done cleaning up' +print('#done cleaning up') cubit.cmd('Display') cubit.cmd('echo on') @@ -88,5 +89,5 @@ # export surface as STL file cubit.cmd('export stl ascii "topo.stl" overwrite') -print '#exporting done' -print '#finished' +print('#exporting done') +print('#finished') diff --git a/EXAMPLES/fault_examples/splay_faults/splay_faults.py b/EXAMPLES/fault_examples/splay_faults/splay_faults.py index 2747047ad..e76c63803 100644 --- a/EXAMPLES/fault_examples/splay_faults/splay_faults.py +++ b/EXAMPLES/fault_examples/splay_faults/splay_faults.py @@ -1,4 +1,5 @@ #!/usr/bin/env python +from __future__ import print_function import cubit import boundary_definition diff --git a/EXAMPLES/fault_examples/tpv102/TPV102.py b/EXAMPLES/fault_examples/tpv102/TPV102.py index fce09305b..7b1fed235 100644 --- a/EXAMPLES/fault_examples/tpv102/TPV102.py +++ b/EXAMPLES/fault_examples/tpv102/TPV102.py @@ -1,6 +1,7 @@ #!python #!/usr/bin/env python # Surendra Nadh Somala, Caltech 2012 +from __future__ import print_function import cubit import cubit2specfem3d diff --git a/EXAMPLES/fault_examples/tpv103/TPV103.py b/EXAMPLES/fault_examples/tpv103/TPV103.py index 3258c0ecc..d17d65f79 100644 --- a/EXAMPLES/fault_examples/tpv103/TPV103.py +++ b/EXAMPLES/fault_examples/tpv103/TPV103.py @@ -1,6 +1,7 @@ #!python #!/usr/bin/env python # Surendra Nadh Somala, Caltech 2012 +from __future__ import print_function import cubit import cubit2specfem3d diff --git a/EXAMPLES/fault_examples/tpv15/TPV15.py b/EXAMPLES/fault_examples/tpv15/TPV15.py index 0d7dd9b1d..3bda58d8b 100644 --- a/EXAMPLES/fault_examples/tpv15/TPV15.py +++ b/EXAMPLES/fault_examples/tpv15/TPV15.py @@ -1,5 +1,6 @@ #!/usr/bin/env python # Percy. Script for TPV14-TPV15 SCEC benchmarks. +from __future__ import print_function import cubit import cubit2specfem3d diff --git a/EXAMPLES/fault_examples/tpv16/TPV16.py b/EXAMPLES/fault_examples/tpv16/TPV16.py index e553ee716..67ceabe79 100644 --- a/EXAMPLES/fault_examples/tpv16/TPV16.py +++ b/EXAMPLES/fault_examples/tpv16/TPV16.py @@ -1,6 +1,7 @@ #!python #!/usr/bin/env python # Surendra Nadh Somala, Caltech 2012 +from __future__ import print_function import cubit import cubit2specfem3d diff --git a/EXAMPLES/fault_examples/tpv5/TPV5.py b/EXAMPLES/fault_examples/tpv5/TPV5.py index 2afa1ebdd..a86b41bde 100755 --- a/EXAMPLES/fault_examples/tpv5/TPV5.py +++ b/EXAMPLES/fault_examples/tpv5/TPV5.py @@ -1,4 +1,6 @@ #!/usr/bin/env python +from __future__ import print_function + import math import os import sys diff --git a/EXAMPLES/homogeneous_halfspace_HEX27_elastic_no_absorbing/block_mesh-anisotropic.py b/EXAMPLES/homogeneous_halfspace_HEX27_elastic_no_absorbing/block_mesh-anisotropic.py index bdd0bcbe3..1f2ca3968 100644 --- a/EXAMPLES/homogeneous_halfspace_HEX27_elastic_no_absorbing/block_mesh-anisotropic.py +++ b/EXAMPLES/homogeneous_halfspace_HEX27_elastic_no_absorbing/block_mesh-anisotropic.py @@ -1,4 +1,5 @@ #!/usr/bin/env python +from __future__ import print_function import cubit cubit.init([""]) diff --git a/EXAMPLES/homogeneous_halfspace_HEX8_elastic_absorbing_Stacey_5sides/create_adjoint_sources.sh b/EXAMPLES/homogeneous_halfspace_HEX8_elastic_absorbing_Stacey_5sides/create_adjoint_sources.sh index ee9179c6a..295e0f66c 100755 --- a/EXAMPLES/homogeneous_halfspace_HEX8_elastic_absorbing_Stacey_5sides/create_adjoint_sources.sh +++ b/EXAMPLES/homogeneous_halfspace_HEX8_elastic_absorbing_Stacey_5sides/create_adjoint_sources.sh @@ -12,6 +12,9 @@ if [ ! -e SEM/xcreate_adjsrc_traveltime ]; then echo "please make xcreate_adjsrc cd SEM/ ./xcreate_adjsrc_traveltime 10.0 25.0 3 X20.DB.BX*.semd +# checks exit code +if [[ $? -ne 0 ]]; then exit 1; fi + if [ ! -e X20.DB.BXZ.semd.adj ]; then echo "error creating adjoint sources, please check..."; exit 1; fi rename .semd.adj .adj *semd.adj diff --git a/EXAMPLES/homogeneous_halfspace_HEX8_elastic_absorbing_Stacey_5sides/create_alpha_kernel_vtk.sh b/EXAMPLES/homogeneous_halfspace_HEX8_elastic_absorbing_Stacey_5sides/create_alpha_kernel_vtk.sh index 5775f1714..ddee428e7 100755 --- a/EXAMPLES/homogeneous_halfspace_HEX8_elastic_absorbing_Stacey_5sides/create_alpha_kernel_vtk.sh +++ b/EXAMPLES/homogeneous_halfspace_HEX8_elastic_absorbing_Stacey_5sides/create_alpha_kernel_vtk.sh @@ -12,6 +12,9 @@ if [ ! -e OUTPUT_FILES/DATABASES_MPI/proc000000_alpha_kernel.bin ]; then echo "p # creates kernel as vtk-file cd bin.xcombine/ ./xcombine_vol_data_vtk 0 3 alpha_kernel ../OUTPUT_FILES/DATABASES_MPI/ ../OUTPUT_FILES/ 1 +# checks exit code +if [[ $? -ne 0 ]]; then exit 1; fi + cd ../ echo diff --git a/EXAMPLES/homogeneous_halfspace_HEX8_elastic_absorbing_Stacey_5sides/create_mesh.py b/EXAMPLES/homogeneous_halfspace_HEX8_elastic_absorbing_Stacey_5sides/create_mesh.py index b29fd6eb9..0b6f9584d 100755 --- a/EXAMPLES/homogeneous_halfspace_HEX8_elastic_absorbing_Stacey_5sides/create_mesh.py +++ b/EXAMPLES/homogeneous_halfspace_HEX8_elastic_absorbing_Stacey_5sides/create_mesh.py @@ -1,9 +1,9 @@ #!/usr/bin/env python - +# # "create_mesh.py" is a script that generates mesh specific to homogenous halfspace example # i.e., a uniform mesh of 134 km x 134 km x 60 km with an element size 3.75 km. # It is not applicable to other examples. - +from __future__ import print_function import os import sys @@ -17,28 +17,28 @@ try: import cubit except ImportError: - print "Error: Importing cubit as python module failed" - print "could not import cubit, please check your PYTHONPATH settings..." - print "" - print "current path: " - print sys.path - print "" - print "try to include path to directory which includes file cubit.py, e.g. /opt/Trelis-15.0/bin/" - print "" + print("Error: Importing cubit as python module failed") + print("could not import cubit, please check your PYTHONPATH settings...") + print("") + print("current path: ") + print(sys.path) + print("") + print("try to include path to directory which includes file cubit.py, e.g. /opt/Trelis-15.0/bin/") + print("") sys.exit("Import cubit failed") -print sys.path +print(sys.path) cubit.init([""]) # gets version string cubit_version = cubit.get_version() -print "version: ",cubit_version +print("version: ",cubit_version) # extracts major number v = cubit_version.split('.') cubit_version_major = int(v[0]) -print "major version number: ",cubit_version_major +print("major version number: ",cubit_version_major) # current work directory cubit.cmd('pwd') @@ -71,9 +71,9 @@ # adds path to geocubit (if not setup yet) sys.path.append('../../CUBIT_GEOCUBIT/') -print "path: " -print sys.path -print "" +print("path: ") +print(sys.path) +print("") # avoids assigning empty blocks cubit.cmd('set duplicate block elements on') @@ -91,17 +91,17 @@ boundary_definition.entities=['face'] boundary_definition.define_bc(boundary_definition.entities,parallel=True) from geocubitlib import cubit2specfem3d - print "" - print "material properties: assigned as block attributes" - print "" + print("") + print("material properties: assigned as block attributes") + print("") # sets the id of the volume block # (volume block starts at id 4) id_block = 1 - print "cubit block:" - print " volume block id = " + str(id_block) - print "" + print("cubit block:") + print(" volume block id = " + str(id_block)) + print("") # Define material properties - print "#### DEFINE MATERIAL PROPERTIES #######################" + print("#### DEFINE MATERIAL PROPERTIES #######################") # elastic material cubit.cmd('block '+str(id_block)+' name "elastic 1" ') # elastic material region cubit.cmd('block '+str(id_block)+' attribute count 7') @@ -119,9 +119,9 @@ #cubit.cmd('block '+str(id_block)+' attribute index 2 1480 ') # vp #cubit.cmd('block '+str(id_block)+' attribute index 3 0 ') # vs #cubit.cmd('block '+str(id_block)+' attribute index 4 1028 ') # rho (ocean salt water density: - print "" - print "exporting to SPECFEM3D-format:" - print "" + print("") + print("exporting to SPECFEM3D-format:") + print("") # Export to SPECFEM3D format cubit2specfem3d.export2SPECFEM3D('MESH/') # backup cubit @@ -129,22 +129,22 @@ cubit.cmd('save as "MESH/meshing.cub" overwrite') else: from geocubitlib import exportlib - print "" - print "exporting to SPECFEM3D-format:" - print "" + print("") + print("exporting to SPECFEM3D-format:") + print("") # Export to SPECFEM3D format # note: exportlib-commands will overwrite material properties exportlib.define_blocks(outdir='MESH/',save_cubfile=True,outfilename='top') exportlib.e2SEM(outdir='MESH/') # Define material properties - print "#### DEFINE MATERIAL PROPERTIES #######################" + print("#### DEFINE MATERIAL PROPERTIES #######################") # elastic material material_cfg=[{'material region':'2','id_block':'1','vp':'2800','vs':'1500','rho':'2300','Qkappa':'9999.0','Qmu':'9999.0','anisotropy_flag':'0'}] # modifies material file nummaterial_velocity_file='MESH/nummaterial_velocity_file' f=open(nummaterial_velocity_file,'w') for block in material_cfg: - print block + print(block) s=block['material region']+' ' s=s+block['id_block']+' ' s=s+block['rho']+' ' diff --git a/EXAMPLES/homogeneous_halfspace_HEX8_elastic_absorbing_Stacey_5sides/make_mesh.sh b/EXAMPLES/homogeneous_halfspace_HEX8_elastic_absorbing_Stacey_5sides/make_mesh.sh index dd6247c9f..1587d3c48 100755 --- a/EXAMPLES/homogeneous_halfspace_HEX8_elastic_absorbing_Stacey_5sides/make_mesh.sh +++ b/EXAMPLES/homogeneous_halfspace_HEX8_elastic_absorbing_Stacey_5sides/make_mesh.sh @@ -12,16 +12,22 @@ echo echo "$geocubit --build_volume --mesh --cfg=homogeneous_halfspace.cfg" echo $geocubit --build_volume --mesh --cfg=homogeneous_halfspace.cfg - +# checks exit code +if [[ $? -ne 0 ]]; then exit 1; fi echo echo "$geocubit --meshfiles=MESH_GEOCUBIT/mesh_vol_0.e --export2SPECFEM3D --SEMoutput=MESH" echo $geocubit --meshfiles=MESH_GEOCUBIT/mesh_vol_0.e --export2SPECFEM3D --SEMoutput=MESH +# checks exit code +if [[ $? -ne 0 ]]; then exit 1; fi echo echo cp -v MESH-default/nummaterial_velocity_file MESH/ +echo +echo "done" +echo diff --git a/EXAMPLES/homogeneous_halfspace_HEX8_elastic_no_absorbing/block_mesh-anisotropic.py b/EXAMPLES/homogeneous_halfspace_HEX8_elastic_no_absorbing/block_mesh-anisotropic.py index d05ee165c..bfcf8d42d 100755 --- a/EXAMPLES/homogeneous_halfspace_HEX8_elastic_no_absorbing/block_mesh-anisotropic.py +++ b/EXAMPLES/homogeneous_halfspace_HEX8_elastic_no_absorbing/block_mesh-anisotropic.py @@ -1,7 +1,12 @@ #!/usr/bin/env python +from __future__ import print_function + +import os +import sys import cubit cubit.init([""]) + try: from geocubitlib import boundary_definition from geocubitlib import cubit2specfem3d @@ -9,10 +14,6 @@ import boundary_definition import cubit2specfem3d - -import os -import sys - # two volumes separating 134000x134000x60000 block horizontally cubit.cmd('reset') cubit.cmd('brick x 67000 y 134000 z 60000') diff --git a/EXAMPLES/layered_halfspace/2lay_mesh_boundary_fig8-nodoubling.py b/EXAMPLES/layered_halfspace/2lay_mesh_boundary_fig8-nodoubling.py index a043f37bc..401777176 100755 --- a/EXAMPLES/layered_halfspace/2lay_mesh_boundary_fig8-nodoubling.py +++ b/EXAMPLES/layered_halfspace/2lay_mesh_boundary_fig8-nodoubling.py @@ -6,6 +6,8 @@ #### specific to the settings of Komatitsch and Tromp 1999, Fig.8 #### Aug 2009 ########################################################################### +from __future__ import print_function + import os import sys @@ -44,7 +46,7 @@ version = cubit.get_version() version_major = int(version.split(".")[0]) version_minor = int(version.split(".")[1]) -print "cubit version: ",version +print("cubit version: ",version) cubit.cmd('reset') cubit.cmd('brick x 134000 y 134000 z 60000') diff --git a/EXAMPLES/layered_halfspace/2lay_mesh_boundary_fig8.py b/EXAMPLES/layered_halfspace/2lay_mesh_boundary_fig8.py index 5734012db..4a63b8453 100755 --- a/EXAMPLES/layered_halfspace/2lay_mesh_boundary_fig8.py +++ b/EXAMPLES/layered_halfspace/2lay_mesh_boundary_fig8.py @@ -6,6 +6,8 @@ #### specific to the settings of Komatitsch and Tromp 1999, Fig.8 #### Aug 2009 ########################################################################### +from __future__ import print_function + import os import sys @@ -34,7 +36,7 @@ version = cubit.get_version() version_major = int(version.split(".")[0]) version_minor = int(version.split(".")[1]) -print "cubit version: ",version +print("cubit version: ",version) cubit.cmd('reset') cubit.cmd('brick x 134000 y 134000 z 60000') diff --git a/EXAMPLES/meshfem3D_examples/cavity/DATA/py/forcesolution_generation.py b/EXAMPLES/meshfem3D_examples/cavity/DATA/py/forcesolution_generation.py index 16b9ab8aa..9f6f0833b 100755 --- a/EXAMPLES/meshfem3D_examples/cavity/DATA/py/forcesolution_generation.py +++ b/EXAMPLES/meshfem3D_examples/cavity/DATA/py/forcesolution_generation.py @@ -1,3 +1,5 @@ +from __future__ import print_function + import os import numpy as np import matplotlib.pyplot as plt diff --git a/EXAMPLES/meshfem3D_examples/cavity/DATA/py/stations_generation.py b/EXAMPLES/meshfem3D_examples/cavity/DATA/py/stations_generation.py index 552894081..3c5f11928 100755 --- a/EXAMPLES/meshfem3D_examples/cavity/DATA/py/stations_generation.py +++ b/EXAMPLES/meshfem3D_examples/cavity/DATA/py/stations_generation.py @@ -1,3 +1,5 @@ +from __future__ import print_function + import os import numpy as np import matplotlib.pyplot as plt diff --git a/EXAMPLES/tomographic_model/tomoblock_mesh.py b/EXAMPLES/tomographic_model/tomoblock_mesh.py index 7e95833ba..9ffcc4a96 100755 --- a/EXAMPLES/tomographic_model/tomoblock_mesh.py +++ b/EXAMPLES/tomographic_model/tomoblock_mesh.py @@ -1,12 +1,12 @@ #!/usr/bin/env python - -import cubit -cubit.init([""]) - +from __future__ import print_function import os import sys +import cubit +cubit.init([""]) + # Creating the volumes cubit.cmd('reset') cubit.cmd('brick x 134000 y 134000 z 60000') @@ -26,9 +26,9 @@ # adds path to geocubit (if not setup yet) sys.path.append('../../CUBIT_GEOCUBIT/') -print "path: " -print sys.path -print "" +print("path: ") +print(sys.path) +print("") try: from geocubitlib import boundary_definition @@ -44,12 +44,12 @@ # sets the id of the volume block # (volume block starts at id 4) id_block = 4 -print "cubit block:" -print " volume block id = " + str(id_block) -print "" +print("cubit block:") +print(" volume block id = " + str(id_block)) +print("") # Define material properties -print "#### DEFINE MATERIAL PROPERTIES #######################" +print("#### DEFINE MATERIAL PROPERTIES #######################") # elastic model cubit.cmd('block '+str(id_block)+' name "elastic tomography_model.xyz 1" ') # elastic material region diff --git a/EXAMPLES/waterlayered_halfspace/make_mesh_waterlayer-default.py b/EXAMPLES/waterlayered_halfspace/make_mesh_waterlayer-default.py index c3e96daca..ef5fdc5b8 100755 --- a/EXAMPLES/waterlayered_halfspace/make_mesh_waterlayer-default.py +++ b/EXAMPLES/waterlayered_halfspace/make_mesh_waterlayer-default.py @@ -8,6 +8,8 @@ # modified version: top layer is water # ########################################################################### +from __future__ import print_function + import os import sys @@ -38,7 +40,7 @@ version = cubit.get_version() version_major = int(version.split(".")[0]) version_minor = int(version.split(".")[1]) -print "cubit version: ",version +print("cubit version: ",version) cubit.cmd('reset') cubit.cmd('brick x 134000 y 134000 z 60000') diff --git a/EXAMPLES/waterlayered_halfspace/make_mesh_waterlayer_verticalboundary.py b/EXAMPLES/waterlayered_halfspace/make_mesh_waterlayer_verticalboundary.py index ba1fa4e07..b8762711c 100755 --- a/EXAMPLES/waterlayered_halfspace/make_mesh_waterlayer_verticalboundary.py +++ b/EXAMPLES/waterlayered_halfspace/make_mesh_waterlayer_verticalboundary.py @@ -1,4 +1,5 @@ #!/usr/bin/env python +from __future__ import print_function SEMoutput='MESH' CUBIToutput='MESH_GEOCUBIT' diff --git a/EXAMPLES/waterlayered_halfspace/waterlayer_mesh_boundary_fig8.py b/EXAMPLES/waterlayered_halfspace/waterlayer_mesh_boundary_fig8.py index 64f28fddc..af8e88bf0 100755 --- a/EXAMPLES/waterlayered_halfspace/waterlayer_mesh_boundary_fig8.py +++ b/EXAMPLES/waterlayered_halfspace/waterlayer_mesh_boundary_fig8.py @@ -1,4 +1,6 @@ #!/usr/bin/env python +from __future__ import print_function + import os import sys diff --git a/EXAMPLES/waterlayered_halfspace/waterlayer_mesh_boundary_vertical.py b/EXAMPLES/waterlayered_halfspace/waterlayer_mesh_boundary_vertical.py index 2e67656c6..0f98a86b3 100755 --- a/EXAMPLES/waterlayered_halfspace/waterlayer_mesh_boundary_vertical.py +++ b/EXAMPLES/waterlayered_halfspace/waterlayer_mesh_boundary_vertical.py @@ -1,4 +1,5 @@ #!/usr/bin/env python +from __future__ import print_function import cubit import boundary_definition diff --git a/EXAMPLES/waterlayered_halfspace/waterlayer_only-nodoubling.py b/EXAMPLES/waterlayered_halfspace/waterlayer_only-nodoubling.py index cb78321bc..a282c3610 100755 --- a/EXAMPLES/waterlayered_halfspace/waterlayer_only-nodoubling.py +++ b/EXAMPLES/waterlayered_halfspace/waterlayer_only-nodoubling.py @@ -1,5 +1,5 @@ #!/usr/bin/env python - +# ########################################################################### #### TNM: This is the mesh generation, adapted from a journal file #### specific to the settings of Komatitsch and Tromp 1999, Fig.8 @@ -9,6 +9,7 @@ # modified version: top layer is water # ########################################################################### +from __future__ import print_function import cubit import boundary_definition diff --git a/src/inverse_problem_for_model/input_output/mesh_tools_mod.f90 b/src/inverse_problem_for_model/input_output/mesh_tools_mod.f90 index e27f8fc50..79ff7af15 100644 --- a/src/inverse_problem_for_model/input_output/mesh_tools_mod.f90 +++ b/src/inverse_problem_for_model/input_output/mesh_tools_mod.f90 @@ -388,7 +388,7 @@ subroutine get_point_in_mesh(x_to_locate, y_to_locate, z_to_locate, iaddx, iaddy ! recompute jacobian for the new point call recompute_jacobian(xelm,yelm,zelm,xi,eta,gamma,x,y,z, & - xixs,xiys,xizs,etaxs,etays,etazs,gammaxs,gammays,gammazs,NGNOD) + xixs,xiys,xizs,etaxs,etays,etazs,gammaxs,gammays,gammazs,NGNOD) ! compute distance to target location dx = - (x - x_target) @@ -418,7 +418,7 @@ subroutine get_point_in_mesh(x_to_locate, y_to_locate, z_to_locate, iaddx, iaddy ! compute final coordinates of point found call recompute_jacobian(xelm,yelm,zelm,xi,eta,gamma,x,y,z, & - xixs,xiys,xizs,etaxs,etays,etazs,gammaxs,gammays,gammazs,NGNOD) + xixs,xiys,xizs,etaxs,etays,etazs,gammaxs,gammays,gammazs,NGNOD) ! store xi,eta,gamma and x,y,z of point found ! note: xi/eta/gamma will be in range [-1,1] diff --git a/utils/Cubit_or_Gmsh/LibGmsh2Specfem_convert_Gmsh_to_Specfem3D.py b/utils/Cubit_or_Gmsh/LibGmsh2Specfem_convert_Gmsh_to_Specfem3D.py index 71bd573bc..f34c5d890 100755 --- a/utils/Cubit_or_Gmsh/LibGmsh2Specfem_convert_Gmsh_to_Specfem3D.py +++ b/utils/Cubit_or_Gmsh/LibGmsh2Specfem_convert_Gmsh_to_Specfem3D.py @@ -22,15 +22,17 @@ # #-------------------------------------------------------------- #-------------------------------------------------------------- - +# # PACKAGES #################################################### +from __future__ import print_function + import sys, string, time from os.path import splitext, isfile try: from numpy import * except ImportError: - print "error: package python-numpy is not installed" + print("error: package python-numpy is not installed") ################################################### @@ -63,8 +65,8 @@ def OuvreGmsh(Dir,Nom): fic=Nom+'.msh' else: - print 'File extension is not correct' - print 'script aborted' + print('File extension is not correct') + print('script aborted') sys.exit() # @@ -127,7 +129,7 @@ def OuvreGmsh(Dir,Nom): Bord_top=Zon ################################################### - print 'Physical Names', PhysCar + print('Physical Names', PhysCar) ################################################### @@ -146,7 +148,7 @@ def OuvreGmsh(Dir,Nom): NbNodes=int(string.split(lignes[PosNodes+1])[0]) # Total number of nodes - print 'Number of nodes: ',NbNodes + print('Number of nodes: ',NbNodes) Nodes=zeros((NbNodes,4),dtype=float) # Array receiving nodes index and coordinates @@ -196,7 +198,7 @@ def OuvreGmsh(Dir,Nom): Ninc2DBordTop, Ninc2DBordBottom, Ninc2DBordxmax, Ninc2DBordxmin, Ninc2DBordymax, Ninc2DBordymin, = 0, 0, 0, 0, 0, 0 - print 'Number of elements: ', NbElements + print('Number of elements: ', NbElements) for Ninc in range(NbElements): @@ -218,8 +220,8 @@ def OuvreGmsh(Dir,Nom): # First case : Surface element - #print 'elem ',Ninc,Pos,TypElem,ZonP,'xmax,xmin,ymax,ymin,top,bottom', \ - # Bord_xmax,Bord_xmin,Bord_ymax,Bord_ymin,Bord_top,Bord_bottom + #print('elem ',Ninc,Pos,TypElem,ZonP,'xmax,xmin,ymax,ymin,top,bottom', \ + # Bord_xmax,Bord_xmin,Bord_ymax,Bord_ymin,Bord_top,Bord_bottom) if TypElem==surfElem: # Get nodes indexes of the surface element @@ -257,12 +259,12 @@ def OuvreGmsh(Dir,Nom): Ninc3D+=1 else: - print "ERROR : wrong element type flag (3 or 5 only)" + print("ERROR : wrong element type flag (3 or 5 only)") # Reduce arrays (exclude zeros elements) - print 'number of elements: xmax,xmin,ymax,ymin,top,bottom = ',Ninc2DBordxmax,Ninc2DBordxmin,Ninc2DBordymax,Ninc2DBordymin, \ - Ninc2DBordTop,Ninc2DBordBottom + print('number of elements: xmax,xmin,ymax,ymin,top,bottom = ', \ + Ninc2DBordxmax,Ninc2DBordxmin,Ninc2DBordymax,Ninc2DBordymin,Ninc2DBordTop,Ninc2DBordBottom) Elements = Elements[:Ninc3D,:] Milieu = Milieu[:Ninc3D,:] @@ -297,10 +299,10 @@ def OuvreGmsh(Dir,Nom): # Nodes in common between 3D current element and boundary rr = set.intersection(nodes3DcurrentElement, NodesBordxmax) if len(rr) != 4: - print "WARNING : wrong 2D boundary element type : ONLY QUADRANGLES" - print "Size of wrong intersection :"+str(len(rr)) - print "Nodes :" - print rr + print("WARNING : wrong 2D boundary element type : ONLY QUADRANGLES") + print("Size of wrong intersection :"+str(len(rr))) + print("Nodes :") + print(rr) sys.exit() else: el = concatenate(([Ct3D+1], list(rr))) @@ -309,10 +311,10 @@ def OuvreGmsh(Dir,Nom): if not set.isdisjoint(nodes3DcurrentElement, NodesBordxmin): rr = set.intersection(nodes3DcurrentElement, NodesBordxmin) if len(rr) != 4: - print "WARNING : wrong 2D boundary element type : ONLY QUADRANGLES" - print "Size of wrong intersection :"+str(len(rr)) - print "Nodes :" - print rr + print("WARNING : wrong 2D boundary element type : ONLY QUADRANGLES") + print("Size of wrong intersection :"+str(len(rr))) + print("Nodes :") + print(rr) sys.exit() else: el = concatenate(([Ct3D+1], list(rr))) @@ -321,10 +323,10 @@ def OuvreGmsh(Dir,Nom): if not set.isdisjoint(nodes3DcurrentElement, NodesBordymax): rr = set.intersection(nodes3DcurrentElement, NodesBordymax) if len(rr) != 4: - print "WARNING : wrong 2D boundary element type : ONLY QUADRANGLES" - print "Size of wrong intersection :"+str(len(rr)) - print "Nodes :" - print rr + print("WARNING : wrong 2D boundary element type : ONLY QUADRANGLES") + print("Size of wrong intersection :"+str(len(rr))) + print("Nodes :") + print(rr) sys.exit() else: el = concatenate(([Ct3D+1], list(rr))) @@ -333,10 +335,10 @@ def OuvreGmsh(Dir,Nom): if not set.isdisjoint(nodes3DcurrentElement, NodesBordymin): rr = set.intersection(nodes3DcurrentElement, NodesBordymin) if len(rr) != 4: - print "WARNING : wrong 2D boundary element type : ONLY QUADRANGLES" - print "Size of wrong intersection :"+str(len(rr)) - print "Nodes :" - print rr + print("WARNING : wrong 2D boundary element type : ONLY QUADRANGLES") + print("Size of wrong intersection :"+str(len(rr))) + print("Nodes :") + print(rr) sys.exit() else: el = concatenate(([Ct3D+1], list(rr))) @@ -345,10 +347,10 @@ def OuvreGmsh(Dir,Nom): if not set.isdisjoint(nodes3DcurrentElement, NodesBordTop): rr = set.intersection(nodes3DcurrentElement, NodesBordTop) if len(rr) != 4: - print "WARNING : wrong 2D boundary element type : ONLY QUADRANGLES" - print "Size of wrong intersection :"+str(len(rr)) - print "Nodes :" - print rr + print("WARNING : wrong 2D boundary element type : ONLY QUADRANGLES") + print("Size of wrong intersection :"+str(len(rr))) + print("Nodes :") + print(rr) sys.exit() else: el = concatenate(([Ct3D+1], list(rr))) @@ -357,10 +359,10 @@ def OuvreGmsh(Dir,Nom): if not set.isdisjoint(nodes3DcurrentElement, NodesBordBottom): rr = set.intersection(nodes3DcurrentElement, NodesBordBottom) if len(rr) != 4: - print "WARNING : wrong 2D boundary element type : ONLY QUADRANGLES" - print "Size of wrong intersection :"+str(len(rr)) - print "Nodes :" - print rr + print("WARNING : wrong 2D boundary element type : ONLY QUADRANGLES") + print("Size of wrong intersection :"+str(len(rr))) + print("Nodes :") + print(rr) sys.exit() else: el = concatenate(([Ct3D+1], list(rr))) diff --git a/utils/Visualization/GMT/grenoble_basin_GMT_movies/mkmpeg4.py b/utils/Visualization/GMT/grenoble_basin_GMT_movies/mkmpeg4.py index 4a0991332..7740e53b1 100755 --- a/utils/Visualization/GMT/grenoble_basin_GMT_movies/mkmpeg4.py +++ b/utils/Visualization/GMT/grenoble_basin_GMT_movies/mkmpeg4.py @@ -1,7 +1,8 @@ #!/usr/bin/env python - +# # Creates an AVI animation from a sequence of images # The input images can be any format other than TIFF with ZLIB compression +from __future__ import print_function import os, sys, re from optparse import OptionParser @@ -59,7 +60,7 @@ def ConvertToSGI(files,verbose): return if verbose: - print 'Converting images to temporary RGB format...' + print('Converting images to temporary RGB format...') for f in files: os.system('convert %s %s.%d.sgi 1> /dev/null 2>&1' % (f, f, pid) ) @@ -78,23 +79,23 @@ def main(): pid = os.getpid() if not out: - print >>sys.stderr,'ERROR: [-o|--out=] output filename required' + print('ERROR: [-o|--out=] output filename required',file=sys.stderr) sys.exit() if not files: - print >>sys.stderr,'ERROR: input files required' + print('ERROR: input files required', file=sys.stderr) sys.exit() if verbose: if divx: - print 'DivX codec selected' + print('DivX codec selected') else: - print 'MS-MPEG4v2 codec selected' + print('MS-MPEG4v2 codec selected') ConvertToSGI(files,verbose) if verbose: - print 'Encoding...pass 1...' + print('Encoding...pass 1...') if divx: opts = "vbitrate=%d:mbd=2:keyint=100:v4mv:vqmin=3:vlelim=-4:vcelim=7:lumi_mask=0.07:dark_mask=0.10:naq:vqcomp=0.7:vqblur=0.2:mpeg_quant" % (bitrate) @@ -106,7 +107,7 @@ def main(): os.system(cmd) if verbose: - print 'Encoding...pass 2...' + print('Encoding...pass 2...') # (Emmanuel Chaljub: this step was not working so I just copy/pasted the encoding pass 1) # cmd = "mencoder -ovc lavc -lavcopts vcodec=mpeg4:vpass=2:%s -mf type=sgi:fps=%d -nosound -o %s mf://\*.%d.sgi" % (opts, fps, out, pid) @@ -127,7 +128,7 @@ def main(): os.system(cmd) if verbose: - print 'Encoding...pass 2...' + print('Encoding...pass 2...') # cmd = "mencoder -ovc lavc -lavcopts vcodec=msmpeg4v2:vpass=2:%s -mf type=sgi:fps=%d -nosound -o %s mf://\*.%d.sgi" % (opts, fps, out, pid) cmd = "mencoder -ovc lavc -lavcopts vcodec=msmpeg4v2:vpass=1:%s -mf type=sgi:fps=%d -nosound -o %s mf://\*.%d.sgi" % (opts, fps, out, pid) @@ -142,7 +143,7 @@ def main(): os.system("rm *.%d.sgi divx2pass.log" % (pid) ) if verbose: - print 'Done' + print('Done') ############################################################## diff --git a/utils/Visualization/Paraview/paraviewpython-example.py b/utils/Visualization/Paraview/paraviewpython-example.py index e8ff110e6..e1fd3155f 100755 --- a/utils/Visualization/Paraview/paraviewpython-example.py +++ b/utils/Visualization/Paraview/paraviewpython-example.py @@ -4,6 +4,7 @@ # usage: ./paraviewpython-example.py alpha_kernel.pvsm [2] # # creates jpg: image*.jpg +from __future__ import print_function import os import sys @@ -20,15 +21,15 @@ filename = str(sys.argv[1]) number = str(sys.argv[2]) else: - print "usage: ./paraviewpython-example.py state-file [counter]" - print " with" - print " state-file - paraview state file, e.g. alpha_kernel.pvsm" - print " counter - (optional) integer counter appended to filename, e.g. 0" + print("usage: ./paraviewpython-example.py state-file [counter]") + print(" with") + print(" state-file - paraview state file, e.g. alpha_kernel.pvsm") + print(" counter - (optional) integer counter appended to filename, e.g. 0") sys.exit(1) #outfile = "paraview_movie." + number outfile = "image" + number -print "file root: ",outfile +print("file root: ",outfile) ## paraview servermanager.Connect() @@ -46,20 +47,20 @@ # sets view size for display #view.ViewSize = [1920,1080] # width x height -print "view size: " + str(view.ViewSize) +print("view size: " + str(view.ViewSize)) ## save as jpeg jpegfilename = outfile + ".jpg" -print "plot to: " + jpegfilename +print("plot to: " + jpegfilename) view.WriteImage(jpegfilename, "vtkJPEGWriter", 1) -print +print("") ## save as png -#print "plot to: " + "image.png" +#print("plot to: " + "image.png") #view.WriteImage("image.png", "vtkPNGWriter",1) ## save as bmp -#print "plot to: " + "image.bmp" +#print("plot to: " + "image.bmp") #view.WriteImage("image.bmp", "vtkBMPWriter",1) diff --git a/utils/Visualization/plot_beachball.py b/utils/Visualization/plot_beachball.py index 8da8b702a..1bb8db1ed 100755 --- a/utils/Visualization/plot_beachball.py +++ b/utils/Visualization/plot_beachball.py @@ -6,6 +6,7 @@ # - obspy : install for example via pip install -U obspy # from __future__ import print_function + import sys import os import datetime diff --git a/utils/combine_result_singals_to_one_hdf5/test_script_convert_signals_to_single_binary.py b/utils/combine_result_singals_to_one_hdf5/test_script_convert_signals_to_single_binary.py index 6cb15a03f..33e0975ce 100644 --- a/utils/combine_result_singals_to_one_hdf5/test_script_convert_signals_to_single_binary.py +++ b/utils/combine_result_singals_to_one_hdf5/test_script_convert_signals_to_single_binary.py @@ -1,3 +1,5 @@ +from __future__ import print_function + import subprocess def main(): diff --git a/utils/combine_result_singals_to_one_hdf5/test_script_open_binarized_file.py b/utils/combine_result_singals_to_one_hdf5/test_script_open_binarized_file.py index 0cf069418..8d73d674c 100644 --- a/utils/combine_result_singals_to_one_hdf5/test_script_open_binarized_file.py +++ b/utils/combine_result_singals_to_one_hdf5/test_script_open_binarized_file.py @@ -1,3 +1,5 @@ +from __future__ import print_function + import h5py def open_bin(fname, id_signal): diff --git a/utils/dynamic_rupture/Visualization/Present_specfem3D_results.py b/utils/dynamic_rupture/Visualization/Present_specfem3D_results.py index 93482d2cc..b5c8389e9 100755 --- a/utils/dynamic_rupture/Visualization/Present_specfem3D_results.py +++ b/utils/dynamic_rupture/Visualization/Present_specfem3D_results.py @@ -1,3 +1,5 @@ +from __future__ import print_function + #from scipy.io import netcdf import numpy as np import sys @@ -41,7 +43,7 @@ ## Files check file_dir = Work_dir + Model_name if not os.path.isdir(file_dir): - print "The directory that contains the data doesn't exist..." + print("The directory that contains the data doesn't exist...") exit() file_list = os.listdir(file_dir) snap_list = [] @@ -52,16 +54,16 @@ time_list.append(int(name.split('Snapshot')[1].split('_')[0])) num_data = len(time_list) time_list = np.asarray(sorted(time_list)) -print -print "The model name is:", Model_name -print "The model has", num_data, "data files." -print "The time list is:", time_list -print +print('') +print("The model name is:", Model_name) +print("The model has", num_data, "data files.") +print("The time list is:", time_list) +print('') if(display_horizontally): - print "Project the fault horizontally." + print("Project the fault horizontally.") else: - print "Project the fault onto a vertical plane (clockwisely rotated ", Rotate, " degree about X-Z plane)." -print + print("Project the fault onto a vertical plane (clockwisely rotated ", Rotate, " degree about X-Z plane).") +print('') #================== # Functions === @@ -110,8 +112,8 @@ def FSEM3D_snapshot(filename): Z = Data[0].Z # if it is a vertical fault (along X or Y axis) if((np.min(X) == np.max(X) or np.min(Y)==np.max(Y)) and display_horizontally): - print "Warning: it is a vertical fault along X or Y axis, the fault will be displayed vertically!!!" - print "Please change display_horizontally to be False and setup the strike of fault to be 0 or 90." + print("Warning: it is a vertical fault along X or Y axis, the fault will be displayed vertically!!!") + print("Please change display_horizontally to be False and setup the strike of fault to be 0 or 90.") exit() if(not display_horizontally): @@ -244,7 +246,7 @@ def FSEM3D_snapshot(filename): plt.savefig("ps/"+Model_name+"-results.pdf",format="pdf") -print "pdf file is saved." +print("pdf file is saved.") ### ### Save data @@ -267,4 +269,4 @@ def FSEM3D_snapshot(filename): output.writelines("\n") output.close() -print "dat file is saved." +print("dat file is saved.") diff --git a/utils/dynamic_rupture/curved_fault_mesh/create_box_mesh.py b/utils/dynamic_rupture/curved_fault_mesh/create_box_mesh.py index 9d535a43e..83e3403c9 100755 --- a/utils/dynamic_rupture/curved_fault_mesh/create_box_mesh.py +++ b/utils/dynamic_rupture/curved_fault_mesh/create_box_mesh.py @@ -1,8 +1,10 @@ #!python # Create 3D mesh files. # Huihui Weng (Geoazur, 2018) - +# # ====================================================================== +from __future__ import print_function + import numpy import os import sys @@ -12,7 +14,7 @@ sys.path.append('/opt/linux64/specfem3d/CUBIT_GEOCUBIT/') import cubit -print "Init CUBIT..." +print("Init CUBIT...") try: # print all the information to the screen. # cubit.init([""]) @@ -94,22 +96,22 @@ # There is no need to change anything below. If you need to change something, # please send me an email. I will try to make it more automatic. # -print "Initial check..." +print("Initial check...") # Initial check if(not os.path.isfile(Int_name) and Interface): - print "The interface data does not exis!!! Please create it in ./Interface." + print("The interface data does not exis!!! Please create it in ./Interface.") exit() elif(os.path.isfile(Int_name) and Interface): - print "Using interface slab: ", Int_name + print("Using interface slab: ", Int_name) else: - print "Using planar fault with strike: ", Strike, " dip: ", Dip, " depth(reference point): ", Dep + print("Using planar fault with strike: ", Strike, " dip: ", Dip, " depth(reference point): ", Dep) if(not os.path.isfile(Top_name) and Topography): - print "The topography data does not exis!!! Please create it in ./Surface." + print("The topography data does not exis!!! Please create it in ./Surface.") elif(os.path.isfile(Top_name) and Topography): - print "Using topography: ", Top_name + print("Using topography: ", Top_name) else: - print "Using planar topography." + print("Using planar topography.") # The name of output mesh file if(Interface and Topography): @@ -131,7 +133,7 @@ output_mesh = output_mesh + "_size" + str(grid_size) + "_" + element_type # Create the journal file for debuging -print "Create journal file..." +print("Create journal file...") j = open(journalFile, 'w') j.write("# Journal file formatting, etc.\n" + \ "# ----------------------------------------------------------------------\n" + \ @@ -255,7 +257,7 @@ j.write("mesh volume {idVol3} {idVol6}\n") j.write("mesh volume all\n") else: - print "Error mesh scheme!" + print("Error mesh scheme!") exit() if(fault_refine_numsplit > 0): j.write("refine surface fault1 NumSplit {0} depth {1}\n".format(fault_refine_numsplit,fault_refine_depth)) @@ -290,7 +292,7 @@ # ================================================== # Read the CUBIT journal and playback it. # ================================================== -print "Playback journal file..." +print("Playback journal file...") with open(journalFile) as f: content = f.readlines() for line in content: @@ -300,8 +302,8 @@ # Save the mesh to txt files # This part is revised from the code of Specfem3D # ================================================== -print "" -print "Convert mesh to Specfem-format..." +print("") +print("Convert mesh to Specfem-format...") os.system('mkdir -p MESH') ## fault surfaces (up/down) @@ -323,15 +325,15 @@ #bottom = [surface_list] #topo = [surface_list] if(len(bottom) == 0 or len(topo) == 0): - print "Fail in obtaining the topo and bottom surfaces." - print "Please change setup topo and bottom surfaces manually." + print("Fail in obtaining the topo and bottom surfaces.") + print("Please change setup topo and bottom surfaces manually.") exit() -print "Xmin surface list: ", xmin -print "Xmax surface list: ", xmax -print "Ymin surface list: ", ymin -print "Ymax surface list: ", ymax -print "Bott surface list: ", bottom -print "Topo surface list: ", topo +print("Xmin surface list: ", xmin) +print("Xmax surface list: ", xmax) +print("Ymin surface list: ", ymin) +print("Ymax surface list: ", ymax) +print("Bott surface list: ", bottom) +print("Topo surface list: ", topo) # define blocks Vol_num = cubit.get_volume_count() @@ -375,7 +377,7 @@ # You need to create fault mesh file in the last, if using hex27. faultA = save_fault_nodes_elements.fault_input(1,Au,Ad) -print "Save created mesh..." +print("Save created mesh...") # Save create directory as given name os.system('rm -rf output/' + output_mesh) os.system('mv MESH output/' + output_mesh) diff --git a/utils/dynamic_rupture/curved_fault_mesh/create_semisphere_mesh.py b/utils/dynamic_rupture/curved_fault_mesh/create_semisphere_mesh.py index 432cf0a4a..6b2939f3f 100755 --- a/utils/dynamic_rupture/curved_fault_mesh/create_semisphere_mesh.py +++ b/utils/dynamic_rupture/curved_fault_mesh/create_semisphere_mesh.py @@ -1,8 +1,10 @@ #!python # Create 3D mesh files. # Huihui Weng (Geoazur, 2018) - +# # ====================================================================== +from __future__ import print_function + import numpy import os import sys @@ -12,7 +14,7 @@ sys.path.append('/opt/linux64/specfem3d/CUBIT_GEOCUBIT/') import cubit -print "Init CUBIT..." +print("Init CUBIT...") try: # print all the information to the screen. cubit.init([""]) @@ -83,22 +85,22 @@ # There is no need to change anything below. If you need to change something, # please send me an email. I will try to make it more automatic. # -print "Initial check..." +print("Initial check...") # Initial check if(not os.path.isfile(Int_name) and Interface): - print "The interface data does not exis!!! Please create it in ./Interface." + print("The interface data does not exis!!! Please create it in ./Interface.") exit() elif(os.path.isfile(Int_name) and Interface): - print "Using interface slab: ", Int_name + print("Using interface slab: ", Int_name) else: - print "Using planar fault with strike: ", Strike, " dip: ", Dip, " depth(reference point): ", Dep + print("Using planar fault with strike: ", Strike, " dip: ", Dip, " depth(reference point): ", Dep) if(not os.path.isfile(Top_name) and Topography): - print "The topography data does not exis!!! Please create it in ./Surface." + print("The topography data does not exis!!! Please create it in ./Surface.") elif(os.path.isfile(Top_name) and Topography): - print "Using topography: ", Top_name + print("Using topography: ", Top_name) else: - print "Using planar topography." + print("Using planar topography.") # The name of output mesh file if(Interface and Topography): @@ -117,7 +119,7 @@ L_cylinder = abs(2 * Lower_cutoff) # Create the journal file for debuging -print "Create journal file..." +print("Create journal file...") j = open(journalFile, 'w') j.write("# Journal file formatting, etc.\n" + \ "# ----------------------------------------------------------------------\n" + \ @@ -220,7 +222,7 @@ j.write("mesh volume {idVol6} \n") j.write("THex Volume all\n") else: - print "Error mesh scheme!" + print("Error mesh scheme!") exit() j.write("# ----------------------------------------------------------------------\n" + \ @@ -251,7 +253,7 @@ # ================================================== # Read the CUBIT journal and playback it. # ================================================== -print "Playback journal file..." +print("Playback journal file...") with open(journalFile) as f: content = f.readlines() for line in content: @@ -261,8 +263,8 @@ # Save the mesh to txt files # This part is revised from the code of Specfem3D # ================================================== -print "" -print "Convert mesh to Specfem-format..." +print("") +print("Convert mesh to Specfem-format...") os.system('mkdir -p MESH') ## fault surfaces (up/down) @@ -281,7 +283,7 @@ if abs(center_point[2]) <= freesur_tolerance: FreesurfID.append(k) -print SpheresurfID,FreesurfID +print(SpheresurfID,FreesurfID) # define blocks Vol_num = cubit.get_volume_count() for i in range(Vol_num): @@ -314,7 +316,7 @@ # You need to create fault mesh file in the last, if using hex27. faultA = save_fault_nodes_elements.fault_input(1,Au,Ad) -print "Save created mesh..." +print("Save created mesh...") # Save create directory as given name os.system('rm -rf output/' + output_mesh) os.system('mv MESH output/' + output_mesh) diff --git a/utils/dynamic_rupture/curved_fault_mesh/python_scripts/Dimension.py b/utils/dynamic_rupture/curved_fault_mesh/python_scripts/Dimension.py index 6ba407336..9e93328b0 100755 --- a/utils/dynamic_rupture/curved_fault_mesh/python_scripts/Dimension.py +++ b/utils/dynamic_rupture/curved_fault_mesh/python_scripts/Dimension.py @@ -1,3 +1,5 @@ +from __future__ import print_function + import sys import numpy @@ -18,5 +20,5 @@ Xnum = Xnum - 2 Ynum = Total_num / Xnum -print Xnum, Ynum +print(Xnum, Ynum) diff --git a/utils/dynamic_rupture/curved_fault_mesh/python_scripts/Smooth_slab.py b/utils/dynamic_rupture/curved_fault_mesh/python_scripts/Smooth_slab.py index 443be0abe..3efef3aff 100755 --- a/utils/dynamic_rupture/curved_fault_mesh/python_scripts/Smooth_slab.py +++ b/utils/dynamic_rupture/curved_fault_mesh/python_scripts/Smooth_slab.py @@ -1,3 +1,5 @@ +from __future__ import print_function + import sys import numpy from scipy.interpolate import griddata @@ -26,6 +28,6 @@ for y in range(Ynum): for x in range(Xnum): if(x%inc == 0 and y%inc == 0 and not numpy.math.isnan(smooth_data[x,y])): - print grid_x[x,y], grid_y[x,y], smooth_data[x,y] + print(grid_x[x,y], grid_y[x,y], smooth_data[x,y]) diff --git a/utils/dynamic_rupture/curved_fault_mesh/python_scripts/playback_jou.py b/utils/dynamic_rupture/curved_fault_mesh/python_scripts/playback_jou.py index 444b8bc01..70fee3210 100644 --- a/utils/dynamic_rupture/curved_fault_mesh/python_scripts/playback_jou.py +++ b/utils/dynamic_rupture/curved_fault_mesh/python_scripts/playback_jou.py @@ -1,3 +1,5 @@ +from __future__ import print_function + import sys # Please set up the path for CUBIT (or Trelis) # Instead, you could set up the path from ~/.bashrc @@ -5,7 +7,7 @@ import cubit -print "Init CUBIT..." +print("Init CUBIT...") try: # output all the information to the screen. #cubit.init([""]) @@ -20,11 +22,11 @@ input_file = P1 -print "Begin to playback the CUBIT script..." +print("Begin to playback the CUBIT script...") ### Read the CUBIT journal and playback it. with open(input_file) as f: content = f.readlines() for line in content: cubit.cmd(line) -print "Finish playback..." +print("Finish playback...") diff --git a/utils/dynamic_rupture/curved_fault_mesh/python_scripts/slab_interface_netsurf.py b/utils/dynamic_rupture/curved_fault_mesh/python_scripts/slab_interface_netsurf.py index dfc0a7da3..5b5f22397 100755 --- a/utils/dynamic_rupture/curved_fault_mesh/python_scripts/slab_interface_netsurf.py +++ b/utils/dynamic_rupture/curved_fault_mesh/python_scripts/slab_interface_netsurf.py @@ -15,8 +15,10 @@ # ---------------------------------------------------------------------- # # PREREQUISITES: numpy - +# # ====================================================================== +from __future__ import print_function + import numpy import sys @@ -28,7 +30,7 @@ P4=sys.argv[i+4] P5=sys.argv[i+5] -print P1,P2,P3,P4,P5 +print(P1,P2,P3,P4,P5) # Define parameters. interfaceFile = P1 numContours = int(P3) @@ -88,7 +90,7 @@ #width = numpy.max(intCoords[:,:,1]) - numpy.min(intCoords[:,:,1]) #offset_x = (numpy.max(intCoords[:,:,0]) + numpy.min(intCoords[:,:,0]))/2.0 #offset_y = (numpy.max(intCoords[:,:,1]) + numpy.min(intCoords[:,:,1]))/2.0 -#print length,width, offset_x, offset_y +#print(length,width, offset_x, offset_y) # # #create surface rectangle width length height width diff --git a/utils/dynamic_rupture/dipping_fault_planar_dip10_kink.py b/utils/dynamic_rupture/dipping_fault_planar_dip10_kink.py index 569e638bc..e3c5a4d83 100644 --- a/utils/dynamic_rupture/dipping_fault_planar_dip10_kink.py +++ b/utils/dynamic_rupture/dipping_fault_planar_dip10_kink.py @@ -2,6 +2,8 @@ # Planar - Dipping Fault. # P. Galvez, ETH-Zurich and J-P Ampuero (Caltech). August 26, 2014. # K. Bai (Caltech). October 2016. Added fault kink (change of dip angle) and re-scaling trick. +from __future__ import print_function + import sys # by STEP3 @@ -24,9 +26,9 @@ def define_bc_topo2(entities,self): topo = self.topo v_list,name_list=define_block() build_block(v_list,name_list) - print entities + print(entities) for entity in entities: - print "##entity: "+str(entity) + print("##entity: "+str(entity)) build_block_side(xmin,entity+'_abs_xmin',obj=entity,id_0=1003) build_block_side(xmax,entity+'_abs_xmax',obj=entity,id_0=1005) build_block_side(ymin,entity+'_abs_ymin',obj=entity,id_0=1004) @@ -62,7 +64,7 @@ def get_global_scale(Xmax,Xmin,Ymax,Ymin,Zmin,h_size): def main(parameter): - print parameter + print(parameter) cubit.init('[]') cubit.cmd('reset') Radian = math.pi/180 @@ -179,7 +181,7 @@ def main(parameter): ##### USER: define material properties ################ cubit.cmd('#### DEFINE MATERIAL PROPERTIES #######################') - print 'we are here' + print('we are here') #for iblock in range(1,11,1): for iblock in range(1,3,1): @@ -193,17 +195,17 @@ def main(parameter): fault = fault_input(1,fu,fd) def usage(): - print "Generate mesh for a planar dipping fault with a kink (change of dip angle)" - print "usage: python dipping_fault_planar_dip10_kink.py L W H A1 A2 ZK DA" - print " L = fault length" - print " W = along-dip width of the fault segment below the kink" - print " H = element size" - print " A1 = shallow dip angle" - print " A2 = deep dip angle" - print " ZK = kink depth" - print " DA = distance from fault edges to absorbing boundaries divided by fault length" - print " All input lengths in km, angles in degree. Output mesh is in m." - print " To change material properties, edit the script in section commented as 'USER'" + print("Generate mesh for a planar dipping fault with a kink (change of dip angle)") + print("usage: python dipping_fault_planar_dip10_kink.py L W H A1 A2 ZK DA") + print(" L = fault length") + print(" W = along-dip width of the fault segment below the kink") + print(" H = element size") + print(" A1 = shallow dip angle") + print(" A2 = deep dip angle") + print(" ZK = kink depth") + print(" DA = distance from fault edges to absorbing boundaries divided by fault length") + print(" All input lengths in km, angles in degree. Output mesh is in m.") + print(" To change material properties, edit the script in section commented as 'USER'") if __name__ == '__main__': if(len(sys.argv) != 8): diff --git a/utils/dynamic_rupture/stepover/change_to_excel_format.py b/utils/dynamic_rupture/stepover/change_to_excel_format.py index 00d876946..ae42f838d 100644 --- a/utils/dynamic_rupture/stepover/change_to_excel_format.py +++ b/utils/dynamic_rupture/stepover/change_to_excel_format.py @@ -1,3 +1,5 @@ +from __future__ import print_function + import sys import os import numpy as np @@ -19,8 +21,5 @@ def convert_to_excel(): f.write('%f, , ,%f,%f,%s,%f\n'%(depth,stress,Vmax,PF,Dmax)) f.close() - - - convert_to_excel() diff --git a/utils/dynamic_rupture/stepover/check_jump_or_not.py b/utils/dynamic_rupture/stepover/check_jump_or_not.py index b4e69e724..900f9bf98 100644 --- a/utils/dynamic_rupture/stepover/check_jump_or_not.py +++ b/utils/dynamic_rupture/stepover/check_jump_or_not.py @@ -1,3 +1,5 @@ +from __future__ import print_function + import numpy as np import sys @@ -11,7 +13,7 @@ def checkjump(time): data = np.fromfile(f1,dtype=np.float32,count=length) tail = np.fromfile(f1,dtype=np.int32,count=1) if(i==8): - print "shear stress",max(data) + print("shear stress",max(data)) with open(prefix+'Snapshot%s_F1.bin'%str(time),'r') as f: for i in range(1,9): @@ -20,12 +22,12 @@ def checkjump(time): data = np.fromfile(f,dtype=np.float32,count=length) tail = np.fromfile(f,dtype=np.int32,count=1) if(i == 2): - print "stepover distance:",max(data) + print("stepover distance:",max(data)) if( i == 4): - print "max displacement fault 1",max(abs(data)) + print("max displacement fault 1",max(abs(data))) Dmax1 = max(abs(data)) if( i == 6): - print "slip veloc on 1",max(abs(data)) + print("slip veloc on 1",max(abs(data))) with open(prefix+'Snapshot%s_F2.bin'%str(time),'r') as f: for i in range(1,9): @@ -34,11 +36,11 @@ def checkjump(time): data = np.fromfile(f,dtype=np.float32,count=length) tail = np.fromfile(f,dtype=np.int32,count=1) if(i == 1): - print "x range",max(data),'-',min(data) + print("x range",max(data),'-',min(data)) if(i == 2): - print "stepover distance:",max(data) + print("stepover distance:",max(data)) if( i == 4 or i == 5): - print "max displacement",max(abs(data)) + print("max displacement",max(abs(data))) if(i == 4): Dmax = max(abs(data)) if(max(abs(data))>2.0): @@ -46,9 +48,9 @@ def checkjump(time): if( i == 6): Vmax = max(abs(data)) if(max(abs(data)) > 0.001): - print "jump happens!",max(abs(data)) + print("jump happens!",max(abs(data))) else: - print "jump fails!",max(abs(data)) + print("jump fails!",max(abs(data))) #this part for checking if the rupture on the first fault goes supershear! with open(prefix+'Snapshot10000_F1.bin','r') as f: for i in range(1,9): @@ -68,7 +70,7 @@ def checkjump(time): return rvalue,Dmax,Vmax,Dmax1,supershear def main(): - print checkjump(20000) + print(checkjump(20000)) if __name__ == "__main__": main() diff --git a/utils/dynamic_rupture/stepover/edit_parfile_fault.py b/utils/dynamic_rupture/stepover/edit_parfile_fault.py index 684bb1585..31c674db5 100644 --- a/utils/dynamic_rupture/stepover/edit_parfile_fault.py +++ b/utils/dynamic_rupture/stepover/edit_parfile_fault.py @@ -1,3 +1,5 @@ +from __future__ import print_function + import sys import check_jump_or_not import re @@ -38,7 +40,7 @@ def setparameter(stress, depth): def run_mpi_job(): p = subprocess.Popen("mpirun -np 90 -perhost 3 -f ~/cpunodes ../bin/xspecfem3D",stdout = subprocess.PIPE,shell = True) - print "mpi job running" + print("mpi job running") p.communicate() def bisect_search(Tmax,Tmin,h): @@ -47,9 +49,9 @@ def bisect_search(Tmax,Tmin,h): if(abs(Tmax - Tmin)<1e5): return (Tmax,Tmin) Tmid = 0.5*(Tmax + Tmin) - print 'now testing T = '+str(Tmid) + print('now testing T = '+str(Tmid)) pf = test(Tmid,h) - print 'test results:'+ str(pf) + print('test results:'+ str(pf)) if(pf): Tmax = Tmid else: @@ -77,7 +79,7 @@ def main(): (a,b) = bisect_search(63e6,69e6,65e3) (a,b) = bisect_search(62e6,68e6,60e3) - print (a,b) + print(a,b) if __name__ == "__main__": diff --git a/utils/dynamic_rupture/stepover/set_uniform_stress.py b/utils/dynamic_rupture/stepover/set_uniform_stress.py index e6332e147..d401e52f1 100644 --- a/utils/dynamic_rupture/stepover/set_uniform_stress.py +++ b/utils/dynamic_rupture/stepover/set_uniform_stress.py @@ -1,4 +1,6 @@ -#!/usr/python +#!/usr/bin/env python +from __future__ import print_function + import os from math import * Mega=1.0e6 diff --git a/utils/dynamic_rupture/subductionzone_mesh/exportmesh.py b/utils/dynamic_rupture/subductionzone_mesh/exportmesh.py index 9d5d15631..7adf2c6fd 100644 --- a/utils/dynamic_rupture/subductionzone_mesh/exportmesh.py +++ b/utils/dynamic_rupture/subductionzone_mesh/exportmesh.py @@ -1,4 +1,6 @@ #!/usr/bin/env python +from __future__ import print_function + import cubit import cubit2specfem3d diff --git a/utils/dynamic_rupture/subductionzone_mesh/process_slab_rotate.py b/utils/dynamic_rupture/subductionzone_mesh/process_slab_rotate.py index e1b22c8c5..68b6a22f1 100644 --- a/utils/dynamic_rupture/subductionzone_mesh/process_slab_rotate.py +++ b/utils/dynamic_rupture/subductionzone_mesh/process_slab_rotate.py @@ -1,6 +1,8 @@ #!/usr/bin/python2.7 #from matplotlib.mlab import griddata #this script generates topology and slab interface in cubit +from __future__ import print_function + import cubit import os import sys