[cig-commits] [commit] devel, master: updated homogeneous_halfspace python scripts for cubit, based on the extra Qkappa column that is needed in nummaterial_velocity_file (1e2c807)

[cig_noreply at geodynamics.org]

Repository : https://github.com/geodynamics/specfem3d

On branches: devel,master
Link       : https://github.com/geodynamics/specfem3d/compare/6026e367984905ab133865f62fa6293b343759b9...47f703851338234f96397e7da9fbff63d8178b8a

>---------------------------------------------------------------

commit 1e2c80752e553649a53e6271b8bf3051f96d2b3c
Author: Carl Tape <carltape at gi.alaska.edu>
Date:   Fri May 17 10:29:20 2013 +0000

    updated homogeneous_halfspace python scripts for cubit, based on the extra Qkappa column that is needed in nummaterial_velocity_file


>---------------------------------------------------------------

1e2c80752e553649a53e6271b8bf3051f96d2b3c
 .../block_mesh.py                                   |  7 ++++---
 .../cubit2specfem3d.py                              | 21 ++++++++++++++-------
 2 files changed, 18 insertions(+), 10 deletions(-)

diff --git a/homogeneous_halfspace_HEX8_elastic_absorbing_Stacey_5sides/block_mesh.py b/homogeneous_halfspace_HEX8_elastic_absorbing_Stacey_5sides/block_mesh.py
index 5730f50..2dcd66a 100755
--- a/homogeneous_halfspace_HEX8_elastic_absorbing_Stacey_5sides/block_mesh.py
+++ b/homogeneous_halfspace_HEX8_elastic_absorbing_Stacey_5sides/block_mesh.py
@@ -31,13 +31,14 @@ cubit.cmd('#### DEFINE MATERIAL PROPERTIES #######################')
 
 
 cubit.cmd('block 1 name "elastic 1" ')        # elastic material region
-cubit.cmd('block 1 attribute count 6')
+cubit.cmd('block 1 attribute count 7')
 cubit.cmd('block 1 attribute index 1 1')      # flag for material: 1 for 1. material
 cubit.cmd('block 1 attribute index 2 2800')   # vp
 cubit.cmd('block 1 attribute index 3 1500')   # vs
 cubit.cmd('block 1 attribute index 4 2300')   # rho
-cubit.cmd('block 1 attribute index 5 9000.0')  # Qmu
-cubit.cmd('block 1 attribute index 6 0 ')      # anisotropy_flag
+cubit.cmd('block 1 attribute index 5 9999.0')  # Qkappa
+cubit.cmd('block 1 attribute index 6 9999.0')  # Qmu
+cubit.cmd('block 1 attribute index 7 0 ')      # anisotropy_flag
 
 
 #cubit.cmd('block 1 name "acoustic 1" ')       # acoustic material region
diff --git a/homogeneous_halfspace_HEX8_elastic_absorbing_Stacey_5sides/cubit2specfem3d.py b/homogeneous_halfspace_HEX8_elastic_absorbing_Stacey_5sides/cubit2specfem3d.py
index 8b29f55..ddc6748 100644
--- a/homogeneous_halfspace_HEX8_elastic_absorbing_Stacey_5sides/cubit2specfem3d.py
+++ b/homogeneous_halfspace_HEX8_elastic_absorbing_Stacey_5sides/cubit2specfem3d.py
@@ -416,9 +416,16 @@ class mesh(object,mesh_tools):
                             if q < 0 :
                               print 'error, q value invalid:', q
                               break
-                            if nattrib == 6:
-                              #anisotropy_flag
-                              ani=cubit.get_block_attribute_value(block,5)
+			    if nattrib >= 6:
+                               #Q_kappa
+                               q2=cubit.get_block_attribute_value(block,5)
+                               # for q to be valid: it must be positive
+                               if q2 < 0 :
+                                  print 'error, q value invalid:', q2
+                                  break
+                               if nattrib == 7:
+                                  #anisotropy_flag
+                                  ani=cubit.get_block_attribute_value(block,6)
                     elif flag < 0:
                         # velocity model
                         vel=name
@@ -431,11 +438,11 @@ class mesh(object,mesh_tools):
                             kind='tomography'
                 else:
                     flag=block
-                    vel,vs,rho,q,ani=(name,0,0,0,0)
+                    vel,vs,rho,q,q2,ani=(name,0,0,0,0,0)
                 block_flag.append(int(flag))
                 block_mat.append(block)
                 if flag > 0:
-                    par=tuple([imaterial,flag,vel,vs,rho,q,ani])
+                    par=tuple([imaterial,flag,vel,vs,rho,q,q2,ani])
                 elif flag < 0:
                     if kind=='interface':
                         par=tuple([imaterial,flag,kind,name,flag_down,flag_up])
@@ -514,8 +521,8 @@ class mesh(object,mesh_tools):
             elif type(vel) != str:
                 # velocity model given as vp,vs,rho,..
                 #format nummaterials file: #material_domain_id #material_id #rho #vp #vs #Q_mu #anisotropy_flag
-                txt='%1i %3i %20f %20f %20f %20f %2i\n' % (properties[0],properties[1],properties[4], \
-                         properties[2],properties[3],properties[5],properties[6])
+                txt='%1i %3i %20f %20f %20f %20f %20f %2i\n' % (properties[0],properties[1],properties[4], \
+                         properties[2],properties[3],properties[5],properties[6],properties[7])
             else:
                 txt='%1i %3i %s \n' % (properties[0],properties[1],properties[2])
         elif flag < 0: