[cig-commits] [commit] devel: fix vol sandwich (4e452ef)
cig_noreply at geodynamics.org
cig_noreply at geodynamics.org
Mon Aug 25 13:29:05 PDT 2014
Repository : https://github.com/geodynamics/specfem3d
On branch : devel
Link : https://github.com/geodynamics/specfem3d/compare/cfaa568dc1b27e21b9743675c5595e76960974cc...17fbe272c6d36b95a968afbb124426a8d10311a6
>---------------------------------------------------------------
commit 4e452ef1495686be147dfde0262a015826fcf752
Author: emanuele casarotti <emanuele.casarotti at ingv.it>
Date: Sat Jul 5 02:48:14 2014 +0200
fix vol sandwich
>---------------------------------------------------------------
4e452ef1495686be147dfde0262a015826fcf752
CUBIT_GEOCUBIT/geocubitlib/volumes.py | 104 ++++++++++++++++++++++++----------
1 file changed, 74 insertions(+), 30 deletions(-)
diff --git a/CUBIT_GEOCUBIT/geocubitlib/volumes.py b/CUBIT_GEOCUBIT/geocubitlib/volumes.py
index 1046213..ca02bee 100644
--- a/CUBIT_GEOCUBIT/geocubitlib/volumes.py
+++ b/CUBIT_GEOCUBIT/geocubitlib/volumes.py
@@ -171,7 +171,10 @@ def layercake_volume_ascii_regulargrid_mpiregularmap(filename=None,verticalsandw
xp=(cfg.xmax-cfg.xmin)/float((nx-1))
yp=(cfg.ymax-cfg.ymin)/float((ny-1))
#
- elev=numpy.zeros([nx,ny,cfg.nz],float)
+ if verticalsandwich:
+ elev=numpy.zeros([nx,ny,cfg.nxvol],float)
+ else:
+ elev=numpy.zeros([nx,ny,cfg.nz],float)
coordx=numpy.zeros([nx,ny],float)
coordy=numpy.zeros([nx,ny],float)
#
@@ -184,8 +187,12 @@ def layercake_volume_ascii_regulargrid_mpiregularmap(filename=None,verticalsandw
ivytot=nelem_chunk_y+1
xstep=xlength #distance between vertex on x
ystep=ylength
- for i in range(0,cfg.nz):
- elev[:,:,i] = cfg.zdepth[i]
+ if verticalsandwich:
+ for i in range(0,cfg.nxvol):
+ elev[:,:,i] = cfg.zdepth[i]
+ else:
+ for i in range(0,cfg.nz):
+ elev[:,:,i] = cfg.xwidth[i]
icoord=0
for iy in range(0,ny):
@@ -219,7 +226,12 @@ def layercake_volume_ascii_regulargrid_mpiregularmap(filename=None,verticalsandw
ivertex=0
#
#create vertex
- for inz in range(0,cfg.nz):
+ if verticalsandwich:
+ nlayer=cfg.nxvol
+ else:
+ nlayer=cfg.nz
+
+ for inz in range(0,nlayer):
if cfg.sea and inz==cfg.nz-1: #sea layer
sealevel=True
bathymetry=False
@@ -233,7 +245,7 @@ def layercake_volume_ascii_regulargrid_mpiregularmap(filename=None,verticalsandw
if cfg.bottomflat and inz == 0: #bottom layer
#
- if cfg.geometry_format == 'ascii':
+ if cfg.geometry_format == 'ascii' and not verticalsandwich:
lv=cubit.get_last_id("vertex")
x_current,y_current=(coordx[nxmin_cpu,nymin_cpu],coordy[nxmin_cpu,nymin_cpu])
@@ -268,38 +280,70 @@ def layercake_volume_ascii_regulargrid_mpiregularmap(filename=None,verticalsandw
cubitcommand= 'create vertex '+ str( x_current )+ ' ' + str( y_current) +' '+ str( cfg.depth_bottom )
cubit.cmd(cubitcommand)
#
+ if verticalsandwich:
+ x_current,y_current=geo2utm(coordx[nxmin_cpu,nymax_cpu],coordy[nxmin_cpu,nymax_cpu],'utm')
+ cubitcommand= 'create vertex '+ str( x_current )+ ' ' + str( y_current) +' '+ str( 0 )
+ cubit.cmd(cubitcommand)
+ #
+ x_current,y_current=geo2utm(coordx[nxmin_cpu,nymin_cpu],coordy[nxmin_cpu,nymin_cpu],'utm')
+ cubitcommand= 'create vertex '+ str( x_current )+ ' ' + str( y_current) +' '+ str( 0 )
+ cubit.cmd(cubitcommand)
+ else:
+ x_current,y_current=geo2utm(coordx[nxmax_cpu,nymax_cpu],coordy[nxmax_cpu,nymax_cpu],'utm')
+ cubitcommand= 'create vertex '+ str( x_current )+ ' ' + str( y_current) +' '+ str( cfg.depth_bottom )
+ cubit.cmd(cubitcommand)
+ #
+ x_current,y_current=geo2utm(coordx[nxmax_cpu,nymin_cpu],coordy[nxmax_cpu,nymin_cpu],'utm')
+ cubitcommand= 'create vertex '+ str( x_current )+ ' ' + str( y_current) +' '+ str( cfg.depth_bottom )
+ cubit.cmd(cubitcommand)
+ #
+ lv2=cubit.get_last_id("vertex")
+ cubitcommand= 'create surface vertex '+str(lv+1)+' to '+str(lv2)
+ cubit.cmd(cubitcommand)
+ #
+ isurf = isurf + 1
+ else:
+ if cfg.geometry_format == 'regmesh':
+ if verticalsandwich:
+ zvertex=cfg.xwidth[inz]
+ lv=cubit.get_last_id("vertex")
+ x_current,y_current=geo2utm(coordx[nxmax_cpu,nymin_cpu],coordy[nxmax_cpu,nymin_cpu],'utm')
+ x_current = x_current + zvertex;
+ cubitcommand= 'create vertex '+ str( x_current )+ ' ' + str( y_current) +' '+ str( cfg.depth_bottom )
+ cubit.cmd(cubitcommand)
+ #
x_current,y_current=geo2utm(coordx[nxmax_cpu,nymax_cpu],coordy[nxmax_cpu,nymax_cpu],'utm')
+ x_current = x_current + zvertex;
cubitcommand= 'create vertex '+ str( x_current )+ ' ' + str( y_current) +' '+ str( cfg.depth_bottom )
cubit.cmd(cubitcommand)
#
+ x_current,y_current=geo2utm(coordx[nxmax_cpu,nymax_cpu],coordy[nxmax_cpu,nymax_cpu],'utm')
+ x_current = x_current + zvertex;
+ cubitcommand= 'create vertex '+ str( x_current )+ ' ' + str( y_current) +' '+ str( 0 )
+ cubit.cmd(cubitcommand)
+ #
x_current,y_current=geo2utm(coordx[nxmax_cpu,nymin_cpu],coordy[nxmax_cpu,nymin_cpu],'utm')
- cubitcommand= 'create vertex '+ str( x_current )+ ' ' + str( y_current) +' '+ str( cfg.depth_bottom )
+ x_current = x_current + zvertex;
+ cubitcommand= 'create vertex '+ str( x_current )+ ' ' + str( y_current) +' '+ str( 0 )
+ cubit.cmd(cubitcommand)
+ else:
+ zvertex=cfg.zdepth[inz]
+ lv=cubit.get_last_id("vertex")
+ x_current,y_current=geo2utm(coordx[nxmin_cpu,nymin_cpu],coordy[nxmin_cpu,nymin_cpu],'utm')
+ cubitcommand= 'create vertex '+ str( x_current )+ ' ' + str( y_current) +' '+ str( zvertex )
cubit.cmd(cubitcommand)
#
- lv2=cubit.get_last_id("vertex")
- cubitcommand= 'create surface vertex '+str(lv+1)+' to '+str(lv2)
+ x_current,y_current=geo2utm(coordx[nxmin_cpu,nymax_cpu],coordy[nxmin_cpu,nymax_cpu],'utm')
+ cubitcommand= 'create vertex '+ str( x_current )+ ' ' + str( y_current) +' '+ str( zvertex )
cubit.cmd(cubitcommand)
#
- isurf = isurf + 1
- else:
- if cfg.geometry_format == 'regmesh':
- zvertex=cfg.zdepth[inz]
- lv=cubit.get_last_id("vertex")
- x_current,y_current=geo2utm(coordx[nxmin_cpu,nymin_cpu],coordy[nxmin_cpu,nymin_cpu],'utm')
- cubitcommand= 'create vertex '+ str( x_current )+ ' ' + str( y_current) +' '+ str( zvertex )
- cubit.cmd(cubitcommand)
- #
- x_current,y_current=geo2utm(coordx[nxmin_cpu,nymax_cpu],coordy[nxmin_cpu,nymax_cpu],'utm')
- cubitcommand= 'create vertex '+ str( x_current )+ ' ' + str( y_current) +' '+ str( zvertex )
- cubit.cmd(cubitcommand)
- #
- x_current,y_current=geo2utm(coordx[nxmax_cpu,nymax_cpu],coordy[nxmax_cpu,nymax_cpu],'utm')
- cubitcommand= 'create vertex '+ str( x_current )+ ' ' + str( y_current) +' '+ str( zvertex )
- cubit.cmd(cubitcommand)
- #
- x_current,y_current=geo2utm(coordx[nxmax_cpu,nymin_cpu],coordy[nxmax_cpu,nymin_cpu],'utm')
- cubitcommand= 'create vertex '+ str( x_current )+ ' ' + str( y_current) +' '+ str( zvertex )
- cubit.cmd(cubitcommand)
+ x_current,y_current=geo2utm(coordx[nxmax_cpu,nymax_cpu],coordy[nxmax_cpu,nymax_cpu],'utm')
+ cubitcommand= 'create vertex '+ str( x_current )+ ' ' + str( y_current) +' '+ str( zvertex )
+ cubit.cmd(cubitcommand)
+ #
+ x_current,y_current=geo2utm(coordx[nxmax_cpu,nymin_cpu],coordy[nxmax_cpu,nymin_cpu],'utm')
+ cubitcommand= 'create vertex '+ str( x_current )+ ' ' + str( y_current) +' '+ str( zvertex )
+ cubit.cmd(cubitcommand)
#
cubitcommand= 'create surface vertex '+str(lv+1)+' '+str(lv+2)+' '+str(lv+3)+' '+str(lv+4)
cubit.cmd(cubitcommand)
@@ -376,10 +420,10 @@ def layercake_volume_ascii_regulargrid_mpiregularmap(filename=None,verticalsandw
#
#!create volume
if not onlysurface:
- if cfg.nz == 1:
+ if nlayer == 1:
nsurface=2
else:
- nsurface=cfg.nz
+ nsurface=nlayer
for inz in range(1,nsurface):
ner=cubit.get_error_count()
create_volume(inz,inz+1,method=cfg.volumecreation_method)
More information about the CIG-COMMITS
mailing list