[cig-commits] r22796 - seismo/3D/SPECFEM3D/trunk/CUBIT
dkomati1 at geodynamics.org
dkomati1 at geodynamics.org
Tue Sep 17 09:14:59 PDT 2013
Author: dkomati1
Date: 2013-09-17 09:14:58 -0700 (Tue, 17 Sep 2013)
New Revision: 22796
Added:
seismo/3D/SPECFEM3D/trunk/CUBIT/LatLongUTMconversion.py
seismo/3D/SPECFEM3D/trunk/CUBIT/__init__.py
seismo/3D/SPECFEM3D/trunk/CUBIT/exportlib.py
seismo/3D/SPECFEM3D/trunk/CUBIT/hex_metric.py
seismo/3D/SPECFEM3D/trunk/CUBIT/local_volume.py
seismo/3D/SPECFEM3D/trunk/CUBIT/menu.py
seismo/3D/SPECFEM3D/trunk/CUBIT/mesh_volume.py
seismo/3D/SPECFEM3D/trunk/CUBIT/quality_log.py
seismo/3D/SPECFEM3D/trunk/CUBIT/read_parameter_cfg.py
seismo/3D/SPECFEM3D/trunk/CUBIT/start.py
seismo/3D/SPECFEM3D/trunk/CUBIT/surfaces.py
seismo/3D/SPECFEM3D/trunk/CUBIT/utilities.py
seismo/3D/SPECFEM3D/trunk/CUBIT/volumes.py
Modified:
seismo/3D/SPECFEM3D/trunk/CUBIT/boundary_definition.py
seismo/3D/SPECFEM3D/trunk/CUBIT/cubit2specfem3d.py
seismo/3D/SPECFEM3D/trunk/CUBIT/run_cubit2specfem3d.py
Log:
committed the new version of cubit2specfem3d.py from GEOCUBIT, and fixed several issues in the calling Python script run_cubit2specfem3d.py
Added: seismo/3D/SPECFEM3D/trunk/CUBIT/LatLongUTMconversion.py
===================================================================
--- seismo/3D/SPECFEM3D/trunk/CUBIT/LatLongUTMconversion.py (rev 0)
+++ seismo/3D/SPECFEM3D/trunk/CUBIT/LatLongUTMconversion.py 2013-09-17 16:14:58 UTC (rev 22796)
@@ -0,0 +1,211 @@
+#!/usr/bin/env python
+
+# Lat Long - UTM, UTM - Lat Long conversions
+
+from math import pi, sin, cos, tan, sqrt
+
+#LatLong- UTM conversion..h
+#definitions for lat/long to UTM and UTM to lat/lng conversions
+#include <string.h>
+
+_deg2rad = pi / 180.0
+_rad2deg = 180.0 / pi
+
+_EquatorialRadius = 2
+_eccentricitySquared = 3
+
+_ellipsoid = [
+# id, Ellipsoid name, Equatorial Radius, square of eccentricity
+# first once is a placeholder only, To allow array indices to match id numbers
+ [ -1, "Placeholder", 0, 0],
+ [ 1, "Airy", 6377563, 0.00667054],
+ [ 2, "Australian National", 6378160, 0.006694542],
+ [ 3, "Bessel 1841", 6377397, 0.006674372],
+ [ 4, "Bessel 1841 (Nambia] ", 6377484, 0.006674372],
+ [ 5, "Clarke 1866", 6378206, 0.006768658],
+ [ 6, "Clarke 1880", 6378249, 0.006803511],
+ [ 7, "Everest", 6377276, 0.006637847],
+ [ 8, "Fischer 1960 (Mercury] ", 6378166, 0.006693422],
+ [ 9, "Fischer 1968", 6378150, 0.006693422],
+ [ 10, "GRS 1967", 6378160, 0.006694605],
+ [ 11, "GRS 1980", 6378137, 0.00669438],
+ [ 12, "Helmert 1906", 6378200, 0.006693422],
+ [ 13, "Hough", 6378270, 0.00672267],
+ [ 14, "International", 6378388, 0.00672267],
+ [ 15, "Krassovsky", 6378245, 0.006693422],
+ [ 16, "Modified Airy", 6377340, 0.00667054],
+ [ 17, "Modified Everest", 6377304, 0.006637847],
+ [ 18, "Modified Fischer 1960", 6378155, 0.006693422],
+ [ 19, "South American 1969", 6378160, 0.006694542],
+ [ 20, "WGS 60", 6378165, 0.006693422],
+ [ 21, "WGS 66", 6378145, 0.006694542],
+ [ 22, "WGS-72", 6378135, 0.006694318],
+ [ 23, "WGS-84", 6378137, 0.00669438]
+]
+
+#Reference ellipsoids derived from Peter H. Dana's website-
+#http://www.utexas.edu/depts/grg/gcraft/notes/datum/elist.html
+#Department of Geography, University of Texas at Austin
+#Internet: pdana at mail.utexas.edu
+#3/22/95
+
+#Source
+#Defense Mapping Agency. 1987b. DMA Technical Report: Supplement to Department of Defense World Geodetic System
+#1984 Technical Report. Part I and II. Washington, DC: Defense Mapping Agency
+
+def LLtoUTM(ReferenceEllipsoid, Lat, Long):
+#converts lat/long to UTM coords. Equations from USGS Bulletin 1532
+#East Longitudes are positive, West longitudes are negative.
+#North latitudes are positive, South latitudes are negative
+#Lat and Long are in decimal degrees
+#Written by Chuck Gantz- chuck.gantz at globalstar.com
+
+ a = _ellipsoid[ReferenceEllipsoid][_EquatorialRadius]
+ eccSquared = _ellipsoid[ReferenceEllipsoid][_eccentricitySquared]
+ k0 = 0.9996
+
+#Make sure the longitude is between -180.00 .. 179.9
+ LongTemp = (Long+180)-int((Long+180)/360)*360-180 # -180.00 .. 179.9
+
+ LatRad = Lat*_deg2rad
+ LongRad = LongTemp*_deg2rad
+
+ ZoneNumber = int((LongTemp + 180)/6) + 1
+
+ if Lat >= 56.0 and Lat < 64.0 and LongTemp >= 3.0 and LongTemp < 12.0:
+ ZoneNumber = 32
+
+ # Special zones for Svalbard
+ if Lat >= 72.0 and Lat < 84.0:
+ if LongTemp >= 0.0 and LongTemp < 9.0:ZoneNumber = 31
+ elif LongTemp >= 9.0 and LongTemp < 21.0: ZoneNumber = 33
+ elif LongTemp >= 21.0 and LongTemp < 33.0: ZoneNumber = 35
+ elif LongTemp >= 33.0 and LongTemp < 42.0: ZoneNumber = 37
+
+ LongOrigin = (ZoneNumber - 1)*6 - 180 + 3 #+3 puts origin in middle of zone
+ LongOriginRad = LongOrigin * _deg2rad
+
+ #compute the UTM Zone from the latitude and longitude
+ UTMZone = "%d%c" % (ZoneNumber, _UTMLetterDesignator(Lat))
+
+ eccPrimeSquared = (eccSquared)/(1-eccSquared)
+ N = a/sqrt(1-eccSquared*sin(LatRad)*sin(LatRad))
+ T = tan(LatRad)*tan(LatRad)
+ C = eccPrimeSquared*cos(LatRad)*cos(LatRad)
+ A = cos(LatRad)*(LongRad-LongOriginRad)
+
+ M = a*((1
+ - eccSquared/4
+ - 3*eccSquared*eccSquared/64
+ - 5*eccSquared*eccSquared*eccSquared/256)*LatRad
+ - (3*eccSquared/8
+ + 3*eccSquared*eccSquared/32
+ + 45*eccSquared*eccSquared*eccSquared/1024)*sin(2*LatRad)
+ + (15*eccSquared*eccSquared/256 + 45*eccSquared*eccSquared*eccSquared/1024)*sin(4*LatRad)
+ - (35*eccSquared*eccSquared*eccSquared/3072)*sin(6*LatRad))
+
+ UTMEasting = (k0*N*(A+(1-T+C)*A*A*A/6
+ + (5-18*T+T*T+72*C-58*eccPrimeSquared)*A*A*A*A*A/120)
+ + 500000.0)
+
+ UTMNorthing = (k0*(M+N*tan(LatRad)*(A*A/2+(5-T+9*C+4*C*C)*A*A*A*A/24
+ + (61
+ -58*T
+ +T*T
+ +600*C
+ -330*eccPrimeSquared)*A*A*A*A*A*A/720)))
+
+ if Lat < 0:
+ UTMNorthing = UTMNorthing + 10000000.0; #10000000 meter offset for southern hemisphere
+ return (UTMZone, UTMEasting, UTMNorthing)
+
+
+def _UTMLetterDesignator(Lat):
+#This routine determines the correct UTM letter designator for the given latitude
+#returns 'Z' if latitude is outside the UTM limits of 84N to 80S
+#Written by Chuck Gantz- chuck.gantz at globalstar.com
+
+ if 84 >= Lat >= 72: return 'X'
+ elif 72 > Lat >= 64: return 'W'
+ elif 64 > Lat >= 56: return 'V'
+ elif 56 > Lat >= 48: return 'U'
+ elif 48 > Lat >= 40: return 'T'
+ elif 40 > Lat >= 32: return 'S'
+ elif 32 > Lat >= 24: return 'R'
+ elif 24 > Lat >= 16: return 'Q'
+ elif 16 > Lat >= 8: return 'P'
+ elif 8 > Lat >= 0: return 'N'
+ elif 0 > Lat >= -8: return 'M'
+ elif -8> Lat >= -16: return 'L'
+ elif -16 > Lat >= -24: return 'K'
+ elif -24 > Lat >= -32: return 'J'
+ elif -32 > Lat >= -40: return 'H'
+ elif -40 > Lat >= -48: return 'G'
+ elif -48 > Lat >= -56: return 'F'
+ elif -56 > Lat >= -64: return 'E'
+ elif -64 > Lat >= -72: return 'D'
+ elif -72 > Lat >= -80: return 'C'
+ else: return 'Z' # if the Latitude is outside the UTM limits
+
+#void UTMtoLL(int ReferenceEllipsoid, const double UTMNorthing, const double UTMEasting, const char* UTMZone,
+# double& Lat, double& Long )
+
+def UTMtoLL(ReferenceEllipsoid, northing, easting, zone):
+
+#converts UTM coords to lat/long. Equations from USGS Bulletin 1532
+#East Longitudes are positive, West longitudes are negative.
+#North latitudes are positive, South latitudes are negative
+#Lat and Long are in decimal degrees.
+#Written by Chuck Gantz- chuck.gantz at globalstar.com
+#Converted to Python by Russ Nelson <nelson at crynwr.com>
+
+ k0 = 0.9996
+ a = _ellipsoid[ReferenceEllipsoid][_EquatorialRadius]
+ eccSquared = _ellipsoid[ReferenceEllipsoid][_eccentricitySquared]
+ e1 = (1-sqrt(1-eccSquared))/(1+sqrt(1-eccSquared))
+ #NorthernHemisphere; //1 for northern hemispher, 0 for southern
+
+ x = easting - 500000.0 #remove 500,000 meter offset for longitude
+ y = northing
+
+ ZoneLetter = zone[-1]
+ ZoneNumber = int(zone[:-1])
+ if ZoneLetter >= 'N':
+ NorthernHemisphere = 1 # point is in northern hemisphere
+ else:
+ NorthernHemisphere = 0 # point is in southern hemisphere
+ y -= 10000000.0 # remove 10,000,000 meter offset used for southern hemisphere
+
+ LongOrigin = (ZoneNumber - 1)*6 - 180 + 3 # +3 puts origin in middle of zone
+
+ eccPrimeSquared = (eccSquared)/(1-eccSquared)
+
+ M = y / k0
+ mu = M/(a*(1-eccSquared/4-3*eccSquared*eccSquared/64-5*eccSquared*eccSquared*eccSquared/256))
+
+ phi1Rad = (mu + (3*e1/2-27*e1*e1*e1/32)*sin(2*mu)
+ + (21*e1*e1/16-55*e1*e1*e1*e1/32)*sin(4*mu)
+ +(151*e1*e1*e1/96)*sin(6*mu))
+ phi1 = phi1Rad*_rad2deg;
+
+ N1 = a/sqrt(1-eccSquared*sin(phi1Rad)*sin(phi1Rad))
+ T1 = tan(phi1Rad)*tan(phi1Rad)
+ C1 = eccPrimeSquared*cos(phi1Rad)*cos(phi1Rad)
+ R1 = a*(1-eccSquared)/pow(1-eccSquared*sin(phi1Rad)*sin(phi1Rad), 1.5)
+ D = x/(N1*k0)
+
+ Lat = phi1Rad - (N1*tan(phi1Rad)/R1)*(D*D/2-(5+3*T1+10*C1-4*C1*C1-9*eccPrimeSquared)*D*D*D*D/24
+ +(61+90*T1+298*C1+45*T1*T1-252*eccPrimeSquared-3*C1*C1)*D*D*D*D*D*D/720)
+ Lat = Lat * _rad2deg
+
+ Long = (D-(1+2*T1+C1)*D*D*D/6+(5-2*C1+28*T1-3*C1*C1+8*eccPrimeSquared+24*T1*T1)
+ *D*D*D*D*D/120)/cos(phi1Rad)
+ Long = LongOrigin + Long * _rad2deg
+ return (Lat, Long)
+
+if __name__ == '__main__':
+ (z, e, n) = LLtoUTM(23, 45.00, -75.00)
+ print z, e, n
+ (lat, lon) = UTMtoLL(23, n, e, z)
+ print lat, lon
+
Added: seismo/3D/SPECFEM3D/trunk/CUBIT/__init__.py
===================================================================
Modified: seismo/3D/SPECFEM3D/trunk/CUBIT/boundary_definition.py
===================================================================
--- seismo/3D/SPECFEM3D/trunk/CUBIT/boundary_definition.py 2013-09-17 00:06:13 UTC (rev 22795)
+++ seismo/3D/SPECFEM3D/trunk/CUBIT/boundary_definition.py 2013-09-17 16:14:58 UTC (rev 22796)
@@ -1,188 +1,103 @@
-#############################################################################
-# boundary_definition.py #
-# this file is part of GEOCUBIT #
-# #
-# Created by Emanuele Casarotti #
-# Copyright (c) 2008 Istituto Nazionale di Geofisica e Vulcanologia #
-# #
-#############################################################################
-# #
-# GEOCUBIT is free software: you can redistribute it and/or modify #
-# it under the terms of the GNU General Public License as published by #
-# the Free Software Foundation, either version 3 of the License, or #
-# (at your option) any later version. #
-# #
-# GEOCUBIT is distributed in the hope that it will be useful, #
-# but WITHOUT ANY WARRANTY; without even the implied warranty of #
-# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the #
-# GNU General Public License for more details. #
-# #
-# You should have received a copy of the GNU General Public License #
-# along with GEOCUBIT. If not, see <http://www.gnu.org/licenses/>. #
-# #
-#############################################################################
-
-def define_absorbing_surf():
- """
- define the absorbing surfaces for a layered topological box where boundary are surfaces parallel to the axis.
- it returns absorbing_surf,absorbing_surf_xmin,absorbing_surf_xmax,absorbing_surf_ymin,absorbing_surf_ymax,absorbing_surf_bottom,topo_surf
- where
- absorbing_surf is the list of all the absorbing boundary surf
- absorbing_surf_xmin is the list of the absorbing boundary surfaces that correnspond to x=xmin
- ...
- absorbing_surf_bottom is the list of the absorbing boundary surfaces that correspond to z=zmin
- """
+try:
+ import start as start
+ cubit = start.start_cubit()
+except:
try:
- cubit.cmd('comment')
+ import cubit
except:
- try:
- import cubit
- cubit.init([""])
- except:
- print 'error importing cubit'
- import sys
- sys.exit()
- absorbing_surf=[]
- absorbing_surf_xmin=[]
- absorbing_surf_xmax=[]
- absorbing_surf_ymin=[]
- absorbing_surf_ymax=[]
- absorbing_surf_bottom=[]
- top_surf=[]
+ print 'error importing cubit, check if cubit is installed'
+ pass
+def list2str(l):
+ if not isinstance(l,list): l=list(l)
+ return ' '.join(str(x) for x in l)
+
+def map_boundary(cpuxmin=0,cpuxmax=1,cpuymin=0,cpuymax=1,cpux=1,cpuy=1):
+ ymin=[]
+ xmax=[]
+ xmin=[]
+ ymax=[]
+ listfull=[]
+ #
+ for ix in range(cpuxmin,cpuxmax):
+ for iy in range(cpuymin,cpuymax):
+ ip=iy*cpux+ix
+ if ix == cpuxmin:
+ xmin.append(ip)
+ if ix == cpuxmax-1:
+ xmax.append(ip)
+ if iy == cpuymin:
+ ymin.append(ip)
+ if iy == cpuymax-1:
+ ymax.append(ip)
+ #
+ listfull.append(ip)
+ return xmin,xmax,ymin,ymax,listfull
+
+
+def define_4side_lateral_surfaces():
list_vol=cubit.parse_cubit_list("volume","all")
- init_n_vol=len(list_vol)
- zmax_box=cubit.get_total_bounding_box("volume",list_vol)[7]
- zmin_box=cubit.get_total_bounding_box("volume",list_vol)[6] #it is the z_min of the box ... box= xmin,xmax,d,ymin,ymax,d,zmin...
- xmin_box=cubit.get_total_bounding_box("volume",list_vol)[0]
- xmax_box=cubit.get_total_bounding_box("volume",list_vol)[1]
- ymin_box=cubit.get_total_bounding_box("volume",list_vol)[3]
- ymax_box=cubit.get_total_bounding_box("volume",list_vol)[4]
-
-
- list_surf=cubit.parse_cubit_list("surface","all")
-# for k in list_surf:
-# center_point = cubit.get_center_point("surface", k)
-# if abs((center_point[0] - xmin_box)/xmin_box) <= 0.005:
-# absorbing_surf_xmin.append(k)
-# absorbing_surf.append(k)
-# elif abs((center_point[0] - xmax_box)/xmax_box) <= 0.005:
-# absorbing_surf_xmax.append(k)
-# absorbing_surf.append(k)
-# elif abs((center_point[1] - ymin_box)/ymin_box) <= 0.005:
-# absorbing_surf_ymin.append(k)
-# absorbing_surf.append(k)
-# elif abs((center_point[1] - ymax_box)/ymax_box) <= 0.005:
-# absorbing_surf_ymax.append(k)
-# absorbing_surf.append(k)
-# elif abs((center_point[2] - zmin_box)/zmin_box) <= 0.005:
-# absorbing_surf_bottom.append(k)
-# absorbing_surf.append(k)
-# else:
-# sbox=cubit.get_bounding_box('surface',k)
-# dz=abs((sbox[7] - zmax_box)/zmax_box)
-# normal=cubit.get_surface_normal(k)
-# zn=normal[2]
-# dn=abs(zn-1)
-# if dz <= 0.001 and dn < 0.2:
-# top_surf.append(k)
+ surf_xmin=[]
+ surf_ymin=[]
+ surf_xmax=[]
+ surf_ymax=[]
+ for id_vol in list_vol:
+ surf_vertical=[]
+ xsurf=[]
+ ysurf=[]
+ tres=0.3
+ lsurf=cubit.get_relatives("volume",id_vol,"surface")
+ for k in lsurf:
+ normal=cubit.get_surface_normal(k)
+ center_point = cubit.get_center_point("surface", k)
+ if normal[2] >= -1*tres and normal[2] <= tres:
+ surf_vertical.append(k)
+ xsurf.append(center_point[0])
+ ysurf.append(center_point[1])
+ surf_xmin.append(surf_vertical[xsurf.index(min(xsurf))])
+ surf_ymin.append(surf_vertical[ysurf.index(min(ysurf))])
+ surf_xmax.append(surf_vertical[xsurf.index(max(xsurf))])
+ surf_ymax.append(surf_vertical[ysurf.index(max(ysurf))])
+ return surf_xmin,surf_ymin,surf_xmax,surf_ymax
- #box lengths
- x_len = abs( xmax_box - xmin_box)
- y_len = abs( ymax_box - ymin_box)
- z_len = abs( zmax_box - zmin_box)
- print '##boundary box: '
- print '## x length: ' + str(x_len)
- print '## y length: ' + str(y_len)
- print '## z length: ' + str(z_len)
- # debug
- print '## xmin: ' + str(xmin_box)
- print '## xmax: ' + str(xmax_box)
- print '## ymin: ' + str(ymin_box)
- print '## ymax: ' + str(ymax_box)
- print '## zmin: ' + str(zmin_box)
- print '## zmax: ' + str(zmax_box)
- ############################################
- ##
- ## tolerance parameters
- ##
- ## modified for surface topography
- ############################################
- absorbing_surface_distance_tolerance=0.1
- topographic_surface_distance_tolerance=0.1
- topographic_surface_normal_tolerance=0.3
+def lateral_boundary_are_absorbing(ip=0,cpuxmin=0,cpuxmax=1,cpuymin=0,cpuymax=1,cpux=1,cpuy=1):
+ #
+ xmin,ymin,xmax,ymax=define_4side_lateral_surfaces()
+ ip_xmin,ip_xmax,ip_ymin,ip_ymax,listfull=map_boundary(cpuxmin,cpuxmax,cpuymin,cpuymax,cpux,cpuy)
+ if not isinstance(ip_xmin, list): ip_xmin=[ip_xmin]
+ if not isinstance(ip_ymin, list): ip_ymin=[ip_ymin]
+ if not isinstance(ip_xmax, list): ip_xmax=[ip_xmax]
+ if not isinstance(ip_ymax, list): ip_ymax=[ip_ymax]
+ #
+ abs_xmin=[]
+ abs_ymin=[]
+ abs_xmax=[]
+ abs_ymax=[]
+ #
+ if ip in ip_xmin:
+ abs_xmin=xmin
+ print 'proc ',ip,' is has absorbing boundary xmin'
+ if ip in ip_ymin:
+ print 'proc ',ip,' is has absorbing boundary ymin'
+ abs_ymin=ymin
+ if ip in ip_xmax:
+ print 'proc ',ip,' is has absorbing boundary xmax'
+ abs_xmax=xmax
+ if ip in ip_ymax:
+ print 'proc ',ip,' is has absorbing boundary ymax'
+ abs_ymax=ymax
+ return abs_xmin,abs_xmax,abs_ymin,abs_ymax
- for k in list_surf:
- center_point = cubit.get_center_point("surface", k)
-
- #debug
- print '##surface: ' + str(k)
- print '## center point: ' + str(center_point)
- if abs((center_point[0] - xmin_box)/x_len) <= absorbing_surface_distance_tolerance:
- #debug
- print '## xmin surface: ' + str(k)
- absorbing_surf_xmin.append(k)
- absorbing_surf.append(k)
- elif abs((center_point[0] - xmax_box)/x_len) <= absorbing_surface_distance_tolerance:
- #debug
- print '## xmax surface: ' + str(k)
- absorbing_surf_xmax.append(k)
- absorbing_surf.append(k)
- elif abs((center_point[1] - ymin_box)/y_len) <= absorbing_surface_distance_tolerance:
- #debug
- print '## ymin surface: ' + str(k)
- absorbing_surf_ymin.append(k)
- absorbing_surf.append(k)
- elif abs((center_point[1] - ymax_box)/y_len) <= absorbing_surface_distance_tolerance:
- #debug
- print '## ymax surface: ' + str(k)
- absorbing_surf_ymax.append(k)
- absorbing_surf.append(k)
- elif abs((center_point[2] - zmin_box)/z_len) <= absorbing_surface_distance_tolerance:
- #debug
- print '## bottom surface: ' + str(k)
- absorbing_surf_bottom.append(k)
- absorbing_surf.append(k)
- else:
- sbox=cubit.get_bounding_box('surface',k)
- dz=abs((sbox[7] - zmax_box)/z_len)
- normal=cubit.get_surface_normal(k)
- zn=normal[2]
- dn=abs(abs(zn)-1)
- #debug
- #print '## surface element: ' + str(k)
- #print '## surface element: zn ' + str(zn)
- #print '## surface element: dn ' + str(dn)
- #print '## surface element: dz ' + str(dz)
- if dz <= topographic_surface_distance_tolerance and dn < topographic_surface_normal_tolerance:
- #debug
- print '## topo surface: ' + str(k)
- top_surf.append(k)
-
- return absorbing_surf,absorbing_surf_xmin,absorbing_surf_xmax,absorbing_surf_ymin,absorbing_surf_ymax,absorbing_surf_bottom,top_surf
-
-def define_absorbing_surf_nopar():
+def define_surf(ip=0,cpuxmin=0,cpuxmax=1,cpuymin=0,cpuymax=1,cpux=1,cpuy=1):
"""
- define the absorbing surfaces for a layered topological box where boundary surfaces are not parallel to the axis.
- it returns absorbing_surf,topo_surf
- where
- absorbing_surf is the list of all the absorbing boundary surf
+ define the surfaces defining the boundaries of the volume
"""
- try:
- cubit.cmd('comment')
- except:
- try:
- import cubit
- cubit.init([""])
- except:
- print 'error importing cubit'
- import sys
- sys.exit()
+ #
from sets import Set
def product(*args, **kwds):
# product('ABCD', 'xy') --> Ax Ay Bx By Cx Cy Dx Dy
@@ -193,16 +108,16 @@
result = [x+[y] for x in result for y in pool]
return result
absorbing_surf=[]
- absorbing_surf_xmin=[]
- absorbing_surf_xmax=[]
- absorbing_surf_ymin=[]
- absorbing_surf_ymax=[]
- absorbing_surf_bottom=[]
+ xmin=[]
+ xmax=[]
+ ymin=[]
+ ymax=[]
+ #
top_surf=[]
+ bottom_surf=[]
list_vol=cubit.parse_cubit_list("volume","all")
- init_n_vol=len(list_vol)
zmax_box=cubit.get_total_bounding_box("volume",list_vol)[7]
- zmin_box=cubit.get_total_bounding_box("volume",list_vol)[6] #it is the z_min of the box ... box= xmin,xmax,d,ymin,ymax,d,zmin...
+ zmin_box=cubit.get_total_bounding_box("volume",list_vol)[6] #it is the z_min of the box ... box= xmin,xmax,d,ymin,ymax,d,zmin...
xmin_box=cubit.get_total_bounding_box("volume",list_vol)[0]
xmax_box=cubit.get_total_bounding_box("volume",list_vol)[1]
ymin_box=cubit.get_total_bounding_box("volume",list_vol)[3]
@@ -210,20 +125,30 @@
list_surf=cubit.parse_cubit_list("surface","all")
lv=[]
for k in list_surf:
- sbox=cubit.get_bounding_box('surface',k)
- dzmax=abs((sbox[7] - zmax_box)/zmax_box)
- dzmin=abs((sbox[6] - zmin_box)/zmin_box)
- normal=cubit.get_surface_normal(k)
- zn=normal[2]
- if dzmax <= 0.001 and zn > 0.7:
- top_surf.append(k)
- list_vertex=cubit.get_relatives('surface',k,'vertex')
- for v in list_vertex:
- valence=cubit.get_valence(v)
- if valence <= 4: #valence 3 is a corner, 4 is a vertex between 2 volumes, > 4 is a vertex not in the boundaries
- lv.append(v)
- elif dzmin <= 0.001 and zn < -0.7:
- absorbing_surf.append(k)
+ sbox=cubit.get_bounding_box('surface',k)
+ if zmax_box == 0 and sbox[7] == 0:
+ dzmax=0
+ elif zmax_box == 0 or sbox[7] == 0:
+ dzmax=abs(sbox[7] - zmax_box)
+ else:
+ dzmax=abs(sbox[7] - zmax_box)/max(abs(sbox[7]),abs(zmax_box))
+ if zmin_box == 0 and sbox[6] == 0:
+ dzmin=0
+ elif zmin_box == 0 or sbox[6] == 0:
+ dzmin=abs(sbox[6] - zmin_box)
+ else:
+ dzmin=abs(sbox[6] - zmin_box)/max(abs(sbox[6]),abs(zmin_box))
+ normal=cubit.get_surface_normal(k)
+ zn=normal[2]
+ if dzmax <= 0.1 and zn > 0.4:
+ top_surf.append(k)
+ list_vertex=cubit.get_relatives('surface',k,'vertex')
+ for v in list_vertex:
+ valence=cubit.get_valence(v)
+ if valence <= 4: #valence 3 is a corner, 4 is a vertex between 2 volumes, > 4 is a vertex not in the boundaries
+ lv.append(v)
+ elif dzmin <= 0.001 and zn < -0.7:
+ bottom_surf.append(k)
lp=[]
combs=product(lv,lv)
for comb in combs:
@@ -233,262 +158,374 @@
if len(c) == 1:
p=cubit.get_center_point("curve",list(c)[0])
lp.append(p)
- for k in list_surf:
+ for k in list_surf:
center_point = cubit.get_center_point("surface", k)
for p in lp:
- if abs((center_point[0] - p[0])/p[0]) <= 0.005 and abs((center_point[1] - p[1])/p[1]) <= 0.005:
- absorbing_surf.append(k)
- break
- return absorbing_surf,top_surf
+ if abs((center_point[0] - p[0])/p[0]) <= 0.001 and abs((center_point[1] - p[1])/p[1]) <= 0.001:
+ absorbing_surf.append(k)
+ break
+ #
+ four_side=True
+ if four_side:
+ xmin,ymin,xmax,ymax=define_4side_lateral_surfaces()
+ abs_xmin,abs_xmax,abs_ymin,abs_ymax=lateral_boundary_are_absorbing(ip,cpuxmin,cpuxmax,cpuymin,cpuymax,cpux,cpuy)
-def define_absorbing_surf_sphere():
- try:
- cubit.cmd('comment')
- except:
- try:
- import cubit
- cubit.init([""])
- except:
- print 'error importing cubit'
- import sys
- sys.exit()
- surf=[]
- list_surf=cubit.parse_cubit_list("surface","all")
- for s in list_surf:
- v=cubit.get_relatives('surface',s,'volume')
- if len(v) == 1:
- surf.append(s)
- return surf
+ return absorbing_surf,abs_xmin,abs_xmax,abs_ymin,abs_ymax,top_surf,bottom_surf,xmin,ymin,xmax,ymax
+
+
def define_block():
- try:
- cubit.cmd('comment')
- except:
- try:
- import cubit
- cubit.init([""])
- except:
- print 'error importing cubit'
- import sys
- sys.exit()
list_vol=cubit.parse_cubit_list("volume","all")
init_n_vol=len(list_vol)
list_name=map(lambda x: 'vol'+x,map(str,list_vol))
return list_vol,list_name
-def build_block(vol_list,name):
+def build_block(vol_list,name,id_0=1):
+ #
from sets import Set
- try:
- cubit.cmd('comment')
- except:
- try:
- import cubit
- cubit.init([""])
- except:
- print 'error importing cubit'
- import sys
- sys.exit()
+ #
block_list=cubit.get_block_id_list()
if len(block_list) > 0:
- id_block=max(block_list)
+ id_block=max(max(block_list),2)+id_0
else:
- id_block=0
+ id_block=2+id_0
for v,n in zip(vol_list,name):
- id_block+=1
- v_other=Set(vol_list)-Set([v])
- #command= 'block '+str(id_block)+' hex in node in vol '+str(v)+' except hex in vol '+str(list(v_other))
- command= 'block '+str(id_block)+' hex in vol '+str(v)
- command = command.replace("["," ").replace("]"," ")
- cubit.cmd(command)
- command = "block "+str(id_block)+" name '"+n+"'"
- cubit.cmd(command)
+ id_block+=1
+ v_other=Set(vol_list)-Set([v])
+ #command= 'block '+str(id_block)+' hex in node in vol '+str(v)+' except hex in vol '+str(list(v_other))
+ command= 'block '+str(id_block)+' hex in vol '+str(v)+' except hex in vol '+str(list(v_other))
+ print command
+ command = command.replace("["," ").replace("]"," ")
+ cubit.cmd(command)
+ command = "block "+str(id_block)+" name '"+n+"'"
+ cubit.cmd(command)
-def build_block_side(surf_list,name,obj='surface'):
- try:
- cubit.cmd('comment')
- except:
- try:
- import cubit
- cubit.init([""])
- except:
- print 'error importing cubit'
- import sys
- sys.exit()
-
+def build_block_side(surf_list,name,obj='surface',id_0=1):
id_nodeset=cubit.get_next_nodeset_id()
- id_block=cubit.get_next_block_id()
-
- # creates command string
+ id_block=id_0
if obj == 'hex':
txt='hex in node in surface'
txt1='block '+str(id_block)+ ' '+ txt +' '+str(list(surf_list))
txt2="block "+str(id_block)+" name '"+name+"'"
txt1=txt1.replace("["," ").replace("]"," ")
+ cubit.cmd(txt1)
+ cubit.cmd(txt2)
elif obj == 'node':
- txt=obj+' in surface'
- txt1= 'nodeset '+str(id_nodeset)+ ' '+ txt +' '+str(list(surf_list))
- txt1 = txt1.replace("["," ").replace("]"," ")
- txt2 = "nodeset "+str(id_nodeset)+" name '"+name+"'"
+ txt=obj+' in surface'
+ txt1= 'nodeset '+str(id_nodeset)+ ' '+ txt +' '+str(list(surf_list))
+ txt1 = txt1.replace("["," ").replace("]"," ")
+ txt2 = "nodeset "+str(id_nodeset)+" name '"+name+"'"
+ cubit.cmd(txt1)
+ cubit.cmd(txt2)
elif obj == 'face' or obj == 'edge':
txt=obj+' in surface'
txt1= 'block '+str(id_block)+ ' '+ txt +' '+str(list(surf_list))
txt1 = txt1.replace("["," ").replace("]"," ")
txt2 = "block "+str(id_block)+" name '"+name+"'"
+ cubit.cmd(txt1)
+ cubit.cmd(txt2)
else:
txt1=''
- # do not execute: block id might be wrong
- print "##block "+str(id_block)+" name '"+name+"_notsupported (only hex,face,edge,node)'"
- txt2=''
+ txt2="block "+str(id_block)+" name "+name+"_notsupported (only hex,face,edge,node)"
+ try:
+ cubit.cmd('comment "'+txt1+'"')
+ cubit.cmd('comment "'+txt2+'"')
+ except:
+ pass
- # executes commands
- print "# command: " + txt1
- print "# command: " + txt2
- cubit.cmd(txt1)
- cubit.cmd(txt2)
-
def define_bc(*args,**keys):
+ id_0=1
+ #
+ #
parallel=keys.get('parallel',True)
closed=keys.get('closed',False)
- if not closed:
- print "##open region"
-
- # model with parallel sides (e.g. a block)
- if parallel:
- surf,xmin,xmax,ymin,ymax,bottom,topo=define_absorbing_surf()
- else:
- # arbitrary geometry
- surf,topo=define_absorbing_surf_nopar()
-
+ ip=keys.get("iproc",0)
+ cpuxmin=keys.get("cpuxmin",0)
+ cpuymin=keys.get("cpuymin",0)
+ cpux=keys.get("cpux",1)
+ cpuy=keys.get("cpuy",1)
+ cpuxmax=keys.get("cpuxmax",cpux)
+ cpuymax=keys.get("cpuymax",cpuy)
+ #
+ if parallel:
+ absorbing_surf,abs_xmin,abs_xmax,abs_ymin,abs_ymax,top_surf,bottom_surf,xmin,ymin,xmax,ymax=define_surf(ip=ip,cpuxmin=cpuxmin,cpuxmax=cpuxmax,cpuymin=cpuymin,cpuymax=cpuymax,cpux=cpux,cpuy=cpuy)
+ #
+ #id_side=cubit.get_next_block_id()
+ #try:
+ # entities=args[0]
+ #except:
+ # entities=['face']
+ #for entity in entities:
+ # build_block_side(top_surf,entity+'_topo',obj=entity,id_0=1) #topo is block 1 so id_0=1
+ # build_block_side(bottom_surf,entity+'_abs_bottom',obj=entity,id_0=2)
+ # build_block_side(xmin,entity+'_xmin',obj=entity,id_0=3)
+ # build_block_side(ymin,entity+'_ymin',obj=entity,id_0=4)
+ # build_block_side(xmax,entity+'_xmax',obj=entity,id_0=5)
+ # build_block_side(ymax,entity+'_ymax',obj=entity,id_0=6)
+ #
+ id_0=cubit.get_next_block_id()
v_list,name_list=define_block()
- build_block(v_list,name_list)
+ build_block(v_list,name_list,id_0)
+ #
+ elif closed:
+ surf=define_absorbing_surf_sphere()
+ v_list,name_list=define_block()
+ build_block(v_list,name_list,id_0)
entities=args[0]
- print entities
+ id_side=1
for entity in entities:
- print "##entity: "+str(entity)
-
- # block for free surface (w/ topography)
- #print '## topo surface block: ' + str(topo)
- if len(topo) == 0:
- print ""
- print "no topo surface found, please create block face_topo manually..."
- print ""
- else:
- build_block_side(topo,entity+'_topo',obj=entity)
-
- # model has parallel sides (e.g. a block model )
- if parallel:
+ build_block_side(surf,entity+'_closedvol',obj=entity,id_0=id_side)
+ id_side=id_side+1
- # blocks for each side
- if len(xmin) == 0:
- print ""
- print "no abs_xmin surface found, please create block manually..."
- print ""
- else:
- build_block_side(xmin,entity+'_abs_xmin',obj=entity)
-
- # blocks for each side
- if len(xmax) == 0:
- print ""
- print "no abs_xmax surface found, please create block manually..."
- print ""
- else:
- build_block_side(xmax,entity+'_abs_xmax',obj=entity)
+#########################################
- # blocks for each side
- if len(ymin) == 0:
- print ""
- print "no abs_xmin surface found, please create block manually..."
- print ""
- else:
- build_block_side(ymin,entity+'_abs_ymin',obj=entity)
- # blocks for each side
- if len(ymax) == 0:
- print ""
- print "no abs_ymax surface found, please create block manually..."
- print ""
- else:
- build_block_side(ymax,entity+'_abs_ymax',obj=entity)
+def extract_bottom_curves(surf_xmin,surf_ymin,surf_xmax,surf_ymax):
+ curve_xmin=[]
+ curve_ymin=[]
+ curve_xmax=[]
+ curve_ymax=[]
+ from sets import Set #UPGRADE.... the sets module is deprecated after python 2.6
+ for s in surf_xmin:
+ lcs=cubit.get_relatives("surface",s,"curve")
+ for lc in lcs:
+ curve_xmin.append(lc)
+ for s in surf_xmax:
+ lcs=cubit.get_relatives("surface",s,"curve")
+ for lc in lcs:
+ curve_xmax.append(lc)
+ for s in surf_ymin:
+ lcs=cubit.get_relatives("surface",s,"curve")
+ for lc in lcs:
+ curve_ymin.append(lc)
+ for s in surf_ymax:
+ lcs=cubit.get_relatives("surface",s,"curve")
+ for lc in lcs:
+ curve_ymax.append(lc)
+ curve_xmin=list(Set(curve_xmin))
+ curve_ymin=list(Set(curve_ymin))
+ curve_xmax=list(Set(curve_xmax))
+ curve_ymax=list(Set(curve_ymax))
+ curve_bottom_xmin=select_bottom_curve(curve_xmin)
+ curve_bottom_ymin=select_bottom_curve(curve_ymin)
+ curve_bottom_xmax=select_bottom_curve(curve_xmax)
+ curve_bottom_ymax=select_bottom_curve(curve_ymax)
+ #
+ return curve_bottom_xmin,curve_bottom_ymin,curve_bottom_xmax,curve_bottom_ymax
- # blocks for each side
- if len(bottom) == 0:
- print ""
- print "no abs_bottom surface found, please create block manually..."
- print ""
- else:
- build_block_side(bottom,entity+'_abs_bottom',obj=entity)
+def select_bottom_curve(lc):
+ z=[]
+ for l in lc:
+ center_point = cubit.get_center_point("curve", l)
+ z.append(center_point[2])
+ result=zip(z,lc)
+ result.sort()
+ return result[0][1]
- # block for all sides together
- # NOTE:
- # this might fail in some CUBIT versions, when elements are already
- # assigned to other blocks
- try:
- build_block_side(surf,entity+'_abs',obj=entity)
- except:
- print "no combined surface with all sides created"
- else:
- # arbitrary geometry
- # puts all elements in single block
- build_block_side(surf,entity+'_abs',obj=entity)
-
+def get_ordered_node_surf(lsurface,icurve):
+ if not isinstance(lsurface,str):
+ lsurf=list2str(lsurface)
+ #
+ if not isinstance(icurve,str):
+ icurvestr=str(icurve)
+ orient_nodes_surf=[]
+ #
+ cubit.cmd('del 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))
+ nnode=len(nodes_ls)
+ #
+ orient=[]
+ cubit.cmd('del group n1')
+ cubit.cmd("group 'n1' add node in curve "+icurvestr)
+ x=cubit.get_bounding_box('curve', icurve)
+ if x[2]>x[5]:
+ idx=0
else:
- print "##closed region"
+ idx=1
+ group1 = cubit.get_id_from_name("n1")
+ nodes1 = list(cubit.get_group_nodes(group1))
+ for n in nodes1:
+ v = cubit.get_nodal_coordinates(n)
+ orient.append(v[idx])
+ result=zip(orient,nodes1)
+ result.sort()
+ nodes2=[c[1] for c in result]
+ for n in nodes2:
+ try:
+ nodes_ls.remove(n)
+ except:
+ pass
+ orient_nodes_surf=orient_nodes_surf+nodes2
+ #
+ while len(orient_nodes_surf) < nnode:
+ cubit.cmd('del group n1')
+ cubit.cmd("group 'n1' add node in edge in node "+str(nodes2).replace('[',' ').replace(']',' '))
+ group1 = cubit.get_id_from_name("n1")
+ nodes1 = list(cubit.get_group_nodes(group1))
+ orient=[]
+ nd=[]
+ for n in nodes1:
+ if n in nodes_ls:
+ v = cubit.get_nodal_coordinates(n)
+ orient.append(v[idx])
+ nd.append(n)
+ result=zip(orient,nd)
+ result.sort()
+ nodes2=[c[1] for c in result]
+ for n in nodes2:
+ try:
+ nodes_ls.remove(n)
+ except:
+ pass
+ orient_nodes_surf=orient_nodes_surf+nodes2
+ #get the vertical curve
+ curve_vertical=[]
+ for s in lsurface:
+ lcs=cubit.get_relatives("surface",s,"curve")
+ for l in lcs:
+ x=cubit.get_bounding_box('curve', l)
+ length=[(x[2],1),(x[5],2),(x[8],3)]
+ length.sort()
+ if length[-1][1] == 3:
+ curve_vertical.append(l)
+ #
+ icurve=list2str(curve_vertical)
+ cubit.cmd('del group curve_vertical')
+ cubit.cmd("group 'curve_vertical' add node in curve "+icurve)
+ group1 = cubit.get_id_from_name('curve_vertical')
+ nodes_curve = list(cubit.get_group_nodes(group1))
+ for n in nodes_curve:
+ try:
+ orient_nodes_surf.remove(n)
+ except:
+ pass
+ #
+ return nodes_curve,orient_nodes_surf
- # model without absorbing boundaries, only one surface, e.g. a sphere
- surf=define_absorbing_surf_sphere()
- v_list,name_list=define_block()
- build_block(v_list,name_list)
- entities=args[0]
- for entity in entities:
- # puts all elements in single block
- build_block_side(surf,entity+'_closedvol',obj=entity)
-
-
-
-
-
-## calling example:
-
-#entities=['face']
-#define_bc(entities,parallel=True)
-#define_bc(entities,parallel=False)
-#define_bc(entities,parallel=False,closed=True)
-
-## block material assigning example:
-
-#block 1 attribute count 5
-#block 2 attribute count 0
-#block 2 name '/prova/interface1'
-#block 2 attribute count 3
-#block 3 attribute count 5
-#block 1 attribute index 2 1500
-#block 1 attribute count 2
-#block 3 attribute index 1 2
-#block 3 attribute index 2 5800
-#block 3 attribute index 3 3900
-#block 3 attribute index 4 1500
-#block 3 attribute index 5 3.5
-#block 1 name 'top'
-#block 3 name 'bottom'
-#block 2 attribute index 1 -1
-#block 2 attribute index 2 2
-#block 2 attribute count 4
-#block 2 attribute index 2 1
-#block 2 attribute index 3 2
-
-# to create block manually:
-#
-# use a commands like:
-#
-# e.g. surface with topography
-# block 2 face in surface 6
-# block 2 name "face_topo"
-#
-# e.g. all surface which are absorbing
-# block 3 face in surface 1 2 3 4 5
-# block 2 name "face_abs"
-
-
+def check_bc(iproc,xmin,xmax,ymin,ymax,cpux,cpuy,cpuxmin,cpuxmax,cpuymin,cpuymax):
+ """
+ boundary=check_bc(iproc,xmin,xmax,ymin,ymax,cpux,cpuy,cpuxmin,cpuxmax,cpuymin,cpuymax)
+ #
+ set the boundary condition during the collecting phase and group the nodes of the vertical surface in groups for the merging phase
+ iproc is the value of the processor
+ xmin,ymin,ymax,ymin are the list of iproc that have at least one absorbing boundary condition
+ """
+ #
+ absorbing_surf,abs_xmin,abs_xmax,abs_ymin,abs_ymax,top_surf,bottom_surf,surf_xmin,surf_ymin,surf_xmax,surf_ymax=define_surf(ip=iproc,cpuxmin=cpuxmin,cpuxmax=cpuxmax,cpuymin=cpuymin,cpuymax=cpuymax,cpux=cpux,cpuy=cpuy)
+ 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_surf,abs_xmin,abs_xmax,abs_ymin,abs_ymax,top_surf,bottom_surf,surf_xmin,surf_ymin,surf_xmax,surf_ymax
+ #
+ #
+ #
+ from sets import Set #UPGRADE.... the sets module is deprecated after python 2.6
+
+ cubit.cmd('set info off')
+ cubit.cmd('set echo off')
+ nodes_curve_ymax,orient_nodes_surf_ymax=get_ordered_node_surf(surf_ymax,curve_bottom_ymax)
+ nodes_curve_xmax,orient_nodes_surf_xmax=get_ordered_node_surf(surf_xmax,curve_bottom_xmax)
+ nodes_curve_ymin,orient_nodes_surf_ymin=get_ordered_node_surf(surf_ymin,curve_bottom_ymin)
+ nodes_curve_xmin,orient_nodes_surf_xmin=get_ordered_node_surf(surf_xmin,curve_bottom_xmin)
+ c_xminymin=Set(nodes_curve_xmin).intersection(nodes_curve_ymin)
+ c_xminymax=Set(nodes_curve_xmin).intersection(nodes_curve_ymax)
+ c_xmaxymin=Set(nodes_curve_xmax).intersection(nodes_curve_ymin)
+ c_xmaxymax=Set(nodes_curve_xmax).intersection(nodes_curve_ymax)
+ cubit.cmd('set info on')
+ cubit.cmd('set echo on')
+
+
+ #
+ orient=[]
+ nd=[]
+ for n in c_xminymin:
+ v = cubit.get_nodal_coordinates(n)
+ orient.append(v[2])
+ nd.append(n)
+ result=zip(orient,nd)
+ result.sort()
+ c_xminymin=[c[1] for c in result]
+ #
+ orient=[]
+ nd=[]
+ for n in c_xminymax:
+ v = cubit.get_nodal_coordinates(n)
+ orient.append(v[2])
+ nd.append(n)
+ result=zip(orient,nd)
+ result.sort()
+ c_xminymax=[c[1] for c in result]
+ #
+ orient=[]
+ nd=[]
+ for n in c_xmaxymin:
+ v = cubit.get_nodal_coordinates(n)
+ orient.append(v[2])
+ nd.append(n)
+ result=zip(orient,nd)
+ result.sort()
+ c_xmaxymin=[c[1] for c in result]
+ #
+ orient=[]
+ nd=[]
+ for n in c_xmaxymax:
+ v = cubit.get_nodal_coordinates(n)
+ orient.append(v[2])
+ nd.append(n)
+ result=zip(orient,nd)
+ result.sort()
+ c_xmaxymax=[c[1] for c in result]
+ #
+ boundary={}
+ boundary['id']=iproc
+ #
+ boundary['nodes_surf_xmin']=orient_nodes_surf_xmin
+ boundary['nodes_surf_xmax']=orient_nodes_surf_xmax
+ boundary['nodes_surf_ymin']=orient_nodes_surf_ymin
+ boundary['nodes_surf_ymax']=orient_nodes_surf_ymax
+ #
+ boundary['node_curve_xminymin']=c_xminymin
+ boundary['node_curve_xminymax']=c_xminymax
+ boundary['node_curve_xmaxymin']=c_xmaxymin
+ boundary['node_curve_xmaxymax']=c_xmaxymax
+ #
+ #
+ #
+ entities=['face']
+ #
+ for entity in entities:
+ if len(abs_xmin) != 0:
+ refname=entity+'_abs_xmin'
+ build_block_side(abs_xmin,refname,obj=entity,id_0=1003)
+ #
+ if len(abs_ymin) != 0:
+ refname=entity+'_abs_ymin'
+ build_block_side(abs_ymin,refname,obj=entity,id_0=1004)
+ #
+ if len(abs_xmax) != 0:
+ refname=entity+'_abs_xmax'
+ build_block_side(abs_xmax,refname,obj=entity,id_0=1005)
+ #
+ if len(abs_ymax) != 0:
+ refname=entity+'_abs_ymax'
+ build_block_side(abs_ymax,refname,obj=entity,id_0=1006)
+ ##
+ refname=entity+'_topo'
+ block=3 #change here..... must be 1 @@@@@@@@@
+ ty=None
+ ty=cubit.get_block_element_type(block)
+ if ty != 'HEX8': cubit.cmd('del block '+str(block))
+ build_block_side(top_surf,refname,obj=entity,id_0=1001)
+ #
+ refname=entity+'_bottom'
+ block=4 #change here..... must be 2 @@@@@@@@
+ ty=None
+ ty=cubit.get_block_element_type(block)
+ if ty != 'HEX8': cubit.cmd('del block '+str(block))
+ build_block_side(bottom_surf,refname,obj=entity,id_0=1002)
+ #
+ #
+ return boundary
\ No newline at end of file
Modified: seismo/3D/SPECFEM3D/trunk/CUBIT/cubit2specfem3d.py
===================================================================
--- seismo/3D/SPECFEM3D/trunk/CUBIT/cubit2specfem3d.py 2013-09-17 00:06:13 UTC (rev 22795)
+++ seismo/3D/SPECFEM3D/trunk/CUBIT/cubit2specfem3d.py 2013-09-17 16:14:58 UTC (rev 22796)
@@ -23,35 +23,32 @@
# #
#############################################################################
#
-#for a complete definition of the format of the mesh in SPECFEM3D_SESAME check the manual (http:/):
+#for a complete definition of the format of the mesh in SPECFEM3D check the manual (http://www.geodynamics.org/cig/software/specfem3d):
#
-#USAGE
+#USAGE
#
#############################################################################
#PREREQUISITE
-#The mesh must be prepared
+#The mesh must be prepared
# automatically using the module boundary_definition (see boundary_definition.py for more information)
-#or
+#or
# manually following the convention:
-# - each material should have a block defined by:
-# material domain_flag (acoustic/elastic/poroelastic)name,flag of the material (integer),p velocity
-# (or the full description: name, flag, vp, vs, rho, Q ... if not present these last 3 parameters will be
-# interpolated by module mat_parameter)
-# - each mesh should have the block definition for the face on the free_surface (topography),
+# - each material should have a block defined by material domain_flag (acoustic/elastic/poroelastic) name,flag of the material (integer),p velocity
+# (or the full description: name, flag, vp, vs, rho, Q ... if not present these last 3 parameters will be interpolated by module mat_parameter)
+# - each mesh should have the block definition for the face on the free_surface (topography),
# the name of this block must be 'face_topo' or you can change the default name in mesh.topo defined in profile.
-# - each mesh should have the block definition for the faces on the absorbing boundaries,
-# one block for each surface with x=Xmin,x=Xmax,y=Ymin,y=Ymax and z=bottom. The names of
-# the blocks should contain the strings "xmin,xmax,ymin,ymax,bottom"
+# - each mesh should have the block definition for the faces on the absorbing boundaries,
+# one block for each surface with x=Xmin,x=Xmax,y=Ymin,y=Ymax and z=bottom. The names of the blocks should contain the strings "xmin,xmax,ymin,ymax,bottom"
#
#############################################################################
#RUN
#In a python script or in the cubit python tab call:
+#
+# export2SPECFEM3D(path_exporting_mesh_SPECFEM3D)
#
-# export2SESAME(path_exporting_mesh_SPECFEM3D_SESAME)
+#the module creates a python class for the mesh: ex. profile=mesh()
+#and it export the files of the mesh needed by the partitioner of SPECFEM3D
#
-#the modele create a python class for the mesh: ex. profile=mesh()
-#and it export the files of the mesh needed by the partitioner of SESAME
-#
#############################################################################
#OUTPUT
#The default output are 11 ASCII files:
@@ -62,61 +59,74 @@
# id_elements id_node1 id_node2 id_node3 id_node4 id_node5 id_node6 id_node7 id_node8
# .....
#
-#__________________________________________________________________________________________
+#__________________________________________________________________________________________
##nodecoord_name='nodes_coords_file' -> the file that contains the coordinates of the nodes of the all mesh
# format:
# number of nodes
# id_node x_coordinate y_coordinate z_coordinate
# .....
#
-#__________________________________________________________________________________________
+#__________________________________________________________________________________________
##material_name='materials_file' -> the file that contains the material flag of the elements
# format:
# id_element flag
# .....
#
-#__________________________________________________________________________________________
+#__________________________________________________________________________________________
##nummaterial_name='nummaterial_velocity_file' -> table of the material properties
# format:
-# flag rho vp vs 0 0 #full definition of the properties, flag > 0
+# #material_domain_id #material_id #rho #vp #vs #Q_mu #anisotropy
# .....
-# flag 'tomography' file_name #for interpolation with tomography
+# #material_domain_id 'tomography' file_name #for interpolation with tomography
# .....
-# flag 'interface' file_name flag_for_the_gll_below_the_interface
-# flag_for_the_gll_above_the_interface #for interpolation with interface
-#__________________________________________________________________________________________
+# #material_domain_id 'interface' file_name flag_for_the_gll_below_the_interface flag_for_the_gll_above_the_interface #for interpolation with interface
+#__________________________________________________________________________________________
##absname='absorbing_surface_file' -> this file contains all the face in all the absorbing boundaries
-##absname_local='absorbing_surface_file'+'_xmin' -> this file contains all the face in the
-# absorbing boundary defined by x=Xmin
-##absname_local='absorbing_surface_file'+'_xmax' -> this file contains all the face in the
-# absorbing boundary defined by x=Xmax
-##absname_local='absorbing_surface_file'+'_ymin' -> this file contains all the face in the
-# absorbing boundary defined by y=Ymin
-##absname_local='absorbing_surface_file'+'_ymax' -> this file contains all the face in the
-# absorbing boundary defined by y=Ymax
-##absname_local='absorbing_surface_file'+'_bottom' -> this file contains all the face in the
-# absorbing boundary defined by z=bottom
+##absname_local='absorbing_surface_file'+'_xmin' -> this file contains all the face in the absorbing boundary defined by x=Xmin
+##absname_local='absorbing_surface_file'+'_xmax' -> this file contains all the face in the absorbing boundary defined by x=Xmax
+##absname_local='absorbing_surface_file'+'_ymin' -> this file contains all the face in the absorbing boundary defined by y=Ymin
+##absname_local='absorbing_surface_file'+'_ymax' -> this file contains all the face in the absorbing boundary defined by y=Ymax
+##absname_local='absorbing_surface_file'+'_bottom' -> this file contains all the face in the absorbing boundary defined by z=bottom
# format:
# number of faces
# id_(element containg the face) id_node1_face id_node2_face id_node3_face id_node4_face
# ....
#
#__________________________________________________________________________________________
-##freename='free_or_absorbing_surface_file_zmax' -> file with the hex on the free surface (usually the topography)
+##freename='free_surface_file' -> file with the hex on the free surface (usually the topography)
# format:
# number of faces
# id_(element containg the face) id_node1_face id_node2_face id_node3_face id_node4_face
#
#__________________________________________________________________________________________
-# it is possible save only one (or more) file singularly: for example if you want only the nodecoord_file
-# call the module mesh.nodescoord_write(full path name)
+#__________________________________________________________________________________________
+##surface='*_surface_file' -> file with the hex on any surface (define by the word 'surface' in the name of the block, ex: moho_surface)
+# optional surfaces, e.g. moho_surface
+# should be created like e.g.:
+# > block 10 face in surface 2
+# > block 10 name 'moho_surface'
#
+#
+# format:
+# number of faces
+# id_(element containg the face) id_node1_face id_node2_face id_node3_face id_node4_face
+#
+#__________________________________________________________________________________________
+# it is possible save only one (or more) file singularly: for example if you want only the nodecoord_file call the module mesh.nodescoord_write(full path name)
+#
#############################################################################
-import cubit
+try:
+ import start as start
+ cubit = start.start_cubit()
+except:
+ try:
+ import cubit
+ except:
+ print 'error importing cubit, check if cubit is installed'
+ pass
class mtools(object):
- """docstring for ciao"""
def __init__(self,frequency,list_surf,list_vp):
super(mtools, self).__init__()
self.frequency = frequency
@@ -128,8 +138,7 @@
def __repr__(self):
txt='Meshing for frequency up to '+str(self.frequency)+'Hz\n'
for surf,vp in zip(self.list_surf,self.list_vp):
- txt=txt+'surface '+str(surf)+', vp ='+str(vp)+' -> size '+str(self.freq2meshsize(vp)[0])\
- +' -> dt '+str(self.freq2meshsize(vp)[0])+'\n'
+ txt=txt+'surface '+str(surf)+', vp ='+str(vp)+' -> size '+str(self.freq2meshsize(vp)[0])+' -> dt '+str(self.freq2meshsize(vp)[0])+'\n'
return txt
def freq2meshsize(self,vp):
velocity=vp*.5
@@ -145,7 +154,7 @@
command = "mesh surf "+str(surf)
cubit.cmd(command)
-class block_tools:
+class block_tools():
def __int__(self):
pass
def create_blocks(self,mesh_entity,list_entity=None,):
@@ -191,16 +200,14 @@
"""Tools for the mesh
#########
dt,edge_dt,freq,edge_freq=seismic_resolution(edges,velocity,bins_d=None,bins_u=None,sidelist=None,ngll=5,np=8)
- Given the velocity of a list of edges, seismic_resolution provides the minimum Dt
- required for the stability condition (and the corrisponding edge).
- Furthermore, given the number of gll point in the element (ngll) and the number
- of GLL point for wavelength, it provide the maximum resolved frequency.
+ Given the velocity of a list of edges, seismic_resolution provides the minimum Dt required for the stability condition (and the corrisponding edge).
+ Furthermore, given the number of gll point in the element (ngll) and the number of GLL point for wavelength, it provide the maximum resolved frequency.
#########
length=edge_length(edge)
return the length of a edge
#########
edge_min,length=edge_min_length(surface)
- given the cubit id of a surface, it return the edge with minimun length
+ given the cubit id of a surface, it return the edge with minimun length
#########
"""
def __int__(self):
@@ -208,10 +215,8 @@
def seismic_resolution(self,edges,velocity,bins_d=None,bins_u=None,sidelist=None):
"""
dt,edge_dt,freq,edge_freq=seismic_resolution(edges,velocity,bins_d=None,bins_u=None,sidelist=None,ngll=5,np=8)
- Given the velocity of a list of edges, seismic_resolution provides the minimum Dt
- required for the stability condition (and the corrisponding edge).
- Furthermore, given the number of gll point in the element (ngll) and the number
- of GLL point for wavelength, it provide the maximum resolved frequency.
+ Given the velocity of a list of edges, seismic_resolution provides the minimum Dt required for the stability condition (and the corrisponding edge).
+ Furthermore, given the number of gll point in the element (ngll) and the number of GLL point for wavelength, it provide the maximum resolved frequency.
"""
ratiostore=1e10
dtstore=1e10
@@ -232,7 +237,6 @@
if ratio >= bin_d and ratio < bin_u:
command = "sideset "+str(side)+" edge "+str(edge)
cubit.cmd(command)
- #print command
break
except:
pass
@@ -251,14 +255,13 @@
def edge_min_length(self,surface):
"""
edge_min,length=edge_min_length(surface)
- given the cubit id of a surface, it return the edge with minimun length
+ given the cubit id of a surface, it return the edge with minimun length
"""
from math import sqrt
self.dmin=99999
edge_store=0
command = "group 'list_edge' add edge in surf "+str(surface)
command = command.replace("["," ").replace("]"," ")
- #print command
cubit.cmd(command)
group=cubit.get_id_from_name("list_edge")
edges=cubit.get_group_edges(group)
@@ -279,8 +282,8 @@
a=[p1[0]-p0[0],p1[1]-p0[1],p1[2]-p0[2]]
b=[p2[0]-p1[0],p2[1]-p1[1],p2[2]-p1[2]]
axb=[a[1]*b[2] - a[2]*b[1], a[2]*b[0] - a[0]*b[2], a[0]*b[1] - a[1]*b[0]]
- dot=0.0
- for i in (0,1,2):
+ dot=0.0
+ for i in (0,1,2):
dot=dot+axb[i]*normal[i]
if dot > 0:
return nodes
@@ -308,7 +311,7 @@
cubit.cmd(command)
command = "sideset "+str(nsideset)+ " name "+ "'ratio-["+str(bin_d)+"_"+str(bin_u)+"['"
cubit.cmd(command)
- nend=cubit.get_next_sideset_id()
+ nend=cubit.get_next_sideset_id()
sidelist=range(nstart,nend)
for block in self.block_mat:
name=cubit.get_exodus_entity_name('block',block)
@@ -320,8 +323,7 @@
es=cubit.get_sub_elements("face", face, 1)
edges=edges+list(es)
edges=Set(edges)
- dtstore,edgedtstore,ratiostore,edgeratiostore=self.seismic_resolution(edges,\
- velocity,bins_d,bins_u,sidelist)
+ dtstore,edgedtstore,ratiostore,edgeratiostore=self.seismic_resolution(edges,velocity,bins_d,bins_u,sidelist)
dt.append(dtstore)
ed_dt.append(edgedtstore)
r.append(ratiostore)
@@ -345,19 +347,22 @@
self.material_name='materials_file'
self.nummaterial_name='nummaterial_velocity_file'
self.absname='absorbing_surface_file'
- self.freename='free_or_absorbing_surface_file_zmax'
+ self.freename='free_surface_file'
self.recname='STATIONS'
- self.face='QUAD4'
- self.face2='SHELL4'
- self.hex='HEX8'
+ version_cubit=float(cubit.get_version())
+ if version_cubit >= 12:
+ self.face='SHELL4'
+ else:
+ self.face='QUAD4'
+ self.hex='HEX'
self.edge='BAR2'
self.topo='face_topo'
self.rec='receivers'
+ self.block_definition()
self.ngll=5
self.percent_gll=0.172
self.point_wavelength=5
- self.block_definition()
- cubit.cmd('compress')
+ cubit.cmd('compress all')
def __repr__(self):
pass
def block_definition(self):
@@ -370,10 +375,10 @@
blocks=cubit.get_block_id_list()
for block in blocks:
name=cubit.get_exodus_entity_name('block',block)
- type=cubit.get_block_element_type(block)
- print block,name,blocks,type,self.hex,self.face
- # block has hexahedral elements (HEX8)
- if type == self.hex:
+ ty=cubit.get_block_element_type(block)
+ #print block,blocks,ty,self.hex,self.face
+ if self.hex in ty:
+ nattrib=cubit.get_block_attribute_count(block)
flag=None
vel=None
vs=None
@@ -381,88 +386,69 @@
q=0
ani=0
# material domain id
- if name.find("acoustic") >= 0 :
+ if "acoustic" in name :
imaterial = 1
- elif name.find("elastic") >= 0 :
+ elif "elastic" in name:
imaterial = 2
- elif name.find("poroelastic") >= 0 :
+ elif "poroelastic" in name:
imaterial = 3
else :
imaterial = 0
- print "block: ",name
- print " could not find appropriate material for this block..."
- print ""
- break
-
- nattrib=cubit.get_block_attribute_count(block)
- if nattrib != 0:
+ #
+ if nattrib > 1:
# material flag:
# positive => material properties,
# negative => interface/tomography domain
flag=int(cubit.get_block_attribute_value(block,0))
if flag > 0 and nattrib >= 2:
- # vp
- vel=cubit.get_block_attribute_value(block,1)
- if nattrib >= 3:
- # vs
- vs=cubit.get_block_attribute_value(block,2)
- if nattrib >= 4:
- #density
- rho=cubit.get_block_attribute_value(block,3)
- if nattrib >= 5:
- #Q_mu
- q=cubit.get_block_attribute_value(block,4)
- # for q to be valid: it must be positive
- if q < 0 :
- print 'error, q value invalid:', q
- break
- if nattrib == 6:
- # only 6 parameters given (skipping Q_kappa ), old format style
- #Q_kappa is 10 times stronger than Q_mu
- q2 = q * 10
- # last entry is anisotropic flag
- ani=cubit.get_block_attribute_value(block,5)
- elif 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)
+ vel=cubit.get_block_attribute_value(block,1)
+ if nattrib >= 3:
+ vs=cubit.get_block_attribute_value(block,2)
+ if nattrib >= 4:
+ rho=cubit.get_block_attribute_value(block,3)
+ if nattrib >= 5:
+ q=cubit.get_block_attribute_value(block,4)
+ # for q to be valid: it must be positive
+ if q < 0 :
+ print 'error, q value invalid:', q
+ break
+ if nattrib == 6:
+ ani=cubit.get_block_attribute_value(block,5)
elif flag < 0:
- # velocity model
vel=name
attrib=cubit.get_block_attribute_value(block,1)
- if attrib == 1:
+ if attrib == 1:
kind='interface'
flag_down=cubit.get_block_attribute_value(block,2)
flag_up=cubit.get_block_attribute_value(block,3)
elif attrib == 2:
kind='tomography'
+ elif nattrib == 1:
+ flag=cubit.get_block_attribute_value(block,0)
+ print 'only 1 attribute ', name,block,flag
+ vel,vs,rho,q,ani=(0,0,0,0,0)
else:
flag=block
- vel,vs,rho,q,q2,ani=(name,0,0,0,0,0)
+ vel,vs,rho,q,ani=(name,0,0,0,0)
block_flag.append(int(flag))
block_mat.append(block)
- if flag > 0:
- par=tuple([imaterial,flag,vel,vs,rho,q,q2,ani])
- elif flag < 0:
+ if flag > 0 and nattrib != 1:
+ par=tuple([imaterial,flag,vel,vs,rho,q,ani])
+ elif flag < 0 and nattrib != 1:
if kind=='interface':
par=tuple([imaterial,flag,kind,name,flag_down,flag_up])
elif kind=='tomography':
par=tuple([imaterial,flag,kind,name])
- elif flag==0:
+ elif flag==0 or nattrib == 1:
par=tuple([imaterial,flag,name])
material[block]=par
- elif (type == self.face) or (type == self.face2) :
- # block has surface elements (QUAD4 or SHELL4)
+ elif ty == self.face: #Stacey condition, we need hex here for pml
block_bc_flag.append(4)
block_bc.append(block)
bc[block]=4 #face has connectivity = 4
- if name == self.topo: topography_face=block
+ if name == self.topo or block == 1001: topography_face=block
+ elif ty == 'SPHERE':
+ pass
else:
# block elements differ from HEX8/QUAD4/SHELL4
print '****************************************'
@@ -473,12 +459,12 @@
print 'please check your block definitions!'
print
print 'only supported types are:'
- print ' HEX8 for volumes'
+ print ' HEX/HEX8 for volumes'
print ' QUAD4 for surface'
print ' SHELL4 for surface'
print '****************************************'
continue
-
+ return None, None,None,None,None,None,None,None
nsets=cubit.get_nodeset_id_list()
if len(nsets) == 0: self.receivers=None
for nset in nsets:
@@ -488,6 +474,14 @@
else:
print 'nodeset '+name+' not defined'
self.receivers=None
+ print block_mat
+ print block_flag
+ print block_bc
+ print block_bc_flag
+ print material
+ print bc
+ print topography_face
+ #
try:
self.block_mat=block_mat
self.block_flag=block_flag
@@ -507,52 +501,61 @@
print bc
print topography
print '****************************************'
- def mat_parameter(self,properties):
- #note: material property acoustic/elastic/poroelastic are defined by the block's name
- print "#material properties:"
+ def mat_parameter(self,properties):
print properties
+ #format nummaterials file: #material_domain_id #material_id #rho #vp #vs #Q_mu #anisotropy_flag
imaterial=properties[0]
flag=properties[1]
if flag > 0:
vel=properties[2]
- if properties[3] is None and type(vel) != str:
- # velocity model scales with given vp value
+ if properties[2] is None and type(vel) != str:
if vel >= 30:
m2km=1000.
else:
m2km=1.
vp=vel/m2km
rho=(1.6612*vp-0.472*vp**2+0.0671*vp**3-0.0043*vp**4+0.000106*vp**4)*m2km
- txt='%1i %3i %20f %20f %20f %1i %1i\n' % (properties[0],properties[1],rho,vel,vel/(3**.5),0,0)
- elif type(vel) != str:
- # velocity model given as vp,vs,rho,..
- #format nummaterials file: #material_domain_id #material_id #rho #vp #vs #Q_mu #anisotropy_flag
- 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])
+ txt='%1i %3i %20f %20f %20f %1i %1i\n' % (properties[0],properties[1],rho,vel,vel/(3**.5),0,0)
+ elif type(vel) != str and vel != 0.:
+ try:
+ q=properties[5]
+ except:
+ q=0.
+ try:
+ ani=properties[6]
+ except:
+ ani=0.
+ #print properties[0],properties[3],properties[1],properties[2],q,ani
+ txt='%1i %3i %20f %20f %20f %20f %20f\n' % (properties[0],properties[1],properties[4],properties[2],properties[3],q,ani)
+ elif type(vel) != str and vel != 0.:
+ helpstring="#material_domain_id #material_id #rho #vp #vs #Q_mu #anisotropy"
+ txt='%1i %3i %s \n' % (properties[0],properties[1],helpstring)
else:
- txt='%1i %3i %s \n' % (properties[0],properties[1],properties[2])
+ helpstring=" --> sintax: #material_domain_id #material_id #rho #vp #vs #Q_mu #anisotropy"
+ txt='%1i %3i %s %s\n' % (properties[0],properties[1],properties[2],helpstring)
elif flag < 0:
if properties[2] == 'tomography':
txt='%1i %3i %s %s\n' % (properties[0],properties[1],properties[2],properties[3])
elif properties[2] == 'interface':
- txt='%1i %3i %s %s %1i %1i\n' % (properties[0],properties[1],properties[2],properties[3],\
- properties[4],properties[5])
+ txt='%1i %3i %s %s %1i %1i\n' % (properties[0],properties[1],properties[2],properties[3],properties[4],properties[5])
+ else:
+ helpstring=" --> sintax: #material_domain_id 'tomography' #file_name "
+ txt='%1i %3i %s %s \n' % (properties[0],properties[1],properties[2],helpstring)
+ #
return txt
def nummaterial_write(self,nummaterial_name):
print 'Writing '+nummaterial_name+'.....'
nummaterial=open(nummaterial_name,'w')
for block in self.block_mat:
#name=cubit.get_exodus_entity_name('block',block)
- nummaterial.write(self.mat_parameter(self.material[block]))
+ nummaterial.write(str(self.mat_parameter(self.material[block])))
nummaterial.close()
- print 'Ok'
def mesh_write(self,mesh_name):
meshfile=open(mesh_name,'w')
print 'Writing '+mesh_name+'.....'
num_elems=cubit.get_hex_count()
print ' number of elements:',str(num_elems)
meshfile.write(str(num_elems)+'\n')
- num_write=0
for block,flag in zip(self.block_mat,self.block_flag):
#print block,flag
hexes=cubit.get_block_hexes(block)
@@ -565,7 +568,6 @@
txt=txt+('%10i %10i %10i %10i %10i %10i %10i %10i\n')% nodes[:]
meshfile.write(txt)
meshfile.close()
- print 'Ok'
def material_write(self,mat_name):
mat=open(mat_name,'w')
print 'Writing '+mat_name+'.....'
@@ -574,7 +576,6 @@
for hexa in hexes:
mat.write(('%10i %10i\n') % (hexa,flag))
mat.close()
- print 'Ok'
def nodescoord_write(self,nodecoord_name):
nodecoord=open(nodecoord_name,'w')
print 'Writing '+nodecoord_name+'.....'
@@ -588,19 +589,17 @@
txt=('%10i %20f %20f %20f\n') % (node,x,y,z)
nodecoord.write(txt)
nodecoord.close()
- print 'Ok'
def free_write(self,freename=None):
- # free surface
cubit.cmd('set info off')
cubit.cmd('set echo off')
cubit.cmd('set journal off')
from sets import Set
normal=(0,0,1)
if not freename: freename=self.freename
- # writes free surface file
+ freehex=open(freename,'w')
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)
@@ -616,17 +615,13 @@
if dic_quads_all.has_key(f):
#print f
nodes=cubit.get_connectivity('face',f)
- if len(nodes) == 0: nodes=cubit.get_connectivity('Face',f)
nodes_ok=self.normal_check(nodes,normal)
- txt='%10i %10i %10i %10i %10i\n' % (h,nodes_ok[0],\
- nodes_ok[1],nodes_ok[2],nodes_ok[3])
+ txt='%10i %10i %10i %10i %10i\n' % (h,nodes_ok[0],nodes_ok[1],nodes_ok[2],nodes_ok[3])
freehex.write(txt)
- freehex.close()
- print 'Ok'
+ freehex.close()
cubit.cmd('set info on')
cubit.cmd('set echo on')
def abs_write(self,absname=None):
- # absorbing boundaries
import re
cubit.cmd('set info off')
cubit.cmd('set echo off')
@@ -634,38 +629,60 @@
from sets import Set
if not absname: absname=self.absname
#
- # loops through all block definitions
+ #
list_hex=cubit.parse_cubit_list('hex','all')
for block,flag in zip(self.block_bc,self.block_bc_flag):
if block != self.topography:
name=cubit.get_exodus_entity_name('block',block)
- print name,block
- absflag=False
+ print ' block name:',name,'id:',block
+ cknormal=True
if re.search('xmin',name):
- filename=absname+'_xmin'
+ print 'xmin'
+ abshex_local=open(absname+'_xmin','w')
normal=(-1,0,0)
elif re.search('xmax',name):
- filename=absname+'_xmax'
+ print "xmax"
+ abshex_local=open(absname+'_xmax','w')
normal=(1,0,0)
elif re.search('ymin',name):
- filename=absname+'_ymin'
+ print "ymin"
+ abshex_local=open(absname+'_ymin','w')
normal=(0,-1,0)
elif re.search('ymax',name):
- filename=absname+'_ymax'
+ print "ymax"
+ abshex_local=open(absname+'_ymax','w')
normal=(0,1,0)
elif re.search('bottom',name):
- filename=absname+'_bottom'
+ print "bottom"
+ abshex_local=open(absname+'_bottom','w')
normal=(0,0,-1)
elif re.search('abs',name):
- #print " ...face_abs - not used so far..."
- filename=absname+'_all'
- absflag=True
+ print "abs all - no implemented yet"
+ cknormal=False
+ abshex_local=open(absname,'w')
else:
- continue
- # opens file
- print 'Writing '+filename+'.....'
- abshex_local=open(filename,'w')
- # gets face elements
+ if block == 1003:
+ print 'xmin'
+ abshex_local=open(absname+'_xmin','w')
+ normal=(-1,0,0)
+ elif block == 1004:
+ print "ymin"
+ abshex_local=open(absname+'_ymin','w')
+ normal=(0,-1,0)
+ elif block == 1005:
+ print "xmax"
+ abshex_local=open(absname+'_xmax','w')
+ normal=(1,0,0)
+ elif block == 1006:
+ print "ymax"
+ abshex_local=open(absname+'_ymax','w')
+ normal=(0,1,0)
+ elif block == 1002:
+ print "bottom"
+ abshex_local=open(absname+'_bottom','w')
+ normal=(0,0,-1)
+ #
+ #
quads_all=cubit.get_block_faces(block)
dic_quads_all=dict(zip(quads_all,quads_all))
print ' number of faces = ',len(quads_all)
@@ -682,19 +699,13 @@
for f in faces:
if dic_quads_all.has_key(f):
nodes=cubit.get_connectivity('face',f)
- if len(nodes) == 0: nodes=cubit.get_connectivity('Face',f)
- if not absflag:
- # checks with specified normal
+ if cknormal:
nodes_ok=self.normal_check(nodes,normal)
- txt='%10i %10i %10i %10i %10i\n' % (h,nodes_ok[0],\
- nodes_ok[1],nodes_ok[2],nodes_ok[3])
else:
- txt='%10i %10i %10i %10i %10i\n' % (h,nodes[0],\
- nodes[1],nodes[2],nodes[3])
+ nodes_ok=nodes
+ txt='%10i %10i %10i %10i %10i\n' % (h,nodes_ok[0],nodes_ok[1],nodes_ok[2],nodes_ok[3])
abshex_local.write(txt)
- # closes file
- abshex_local.close()
- print 'Ok'
+ abshex_local.close()
cubit.cmd('set info on')
cubit.cmd('set echo on')
def surface_write(self,pathdir=None):
@@ -735,13 +746,11 @@
for f in faces:
if dic_quads_all.has_key(f):
nodes=cubit.get_connectivity('face',f)
- if len(nodes) == 0: nodes=cubit.get_connectivity('Face',f)
txt='%10i %10i %10i %10i %10i\n' % (h,nodes[0],\
nodes[1],nodes[2],nodes[3])
surfhex_local.write(txt)
# closes file
surfhex_local.close()
- print 'Ok'
def rec_write(self,recname):
print 'Writing '+self.recname+'.....'
recfile=open(self.recname,'w')
@@ -750,43 +759,36 @@
x,y,z=cubit.get_nodal_coordinates(n)
recfile.write('ST%i XX %20f %20f 0.0 0.0 \n' % (i,x,z))
recfile.close()
- print 'Ok'
def write(self,path=''):
cubit.cmd('set info off')
cubit.cmd('set echo off')
cubit.cmd('set journal off')
+ cubit.cmd('compress all')
if len(path) != 0:
if path[-1] != '/': path=path+'/'
- # mesh file
self.mesh_write(path+self.mesh_name)
- # mesh material
self.material_write(path+self.material_name)
- # mesh coordinates
self.nodescoord_write(path+self.nodecoord_name)
- # material definitions
- self.nummaterial_write(path+self.nummaterial_name)
- # free surface: face_top
self.free_write(path+self.freename)
- # absorbing surfaces: abs_***
self.abs_write(path+self.absname)
+ self.nummaterial_write(path+self.nummaterial_name)
# any other surfaces: ***surface***
self.surface_write(path)
- # receivers
if self.receivers: self.rec_write(path+self.recname)
cubit.cmd('set info on')
cubit.cmd('set echo on')
-def export2SESAME(path_exporting_mesh_SPECFEM3D_SESAME):
- cubit.cmd('set info on')
- cubit.cmd('set echo on')
+def export2SPECFEM3D(path_exporting_mesh_SPECFEM3D='.'):
sem_mesh=mesh()
- sem_mesh.write(path=path_exporting_mesh_SPECFEM3D_SESAME)
+ #sem_mesh.block_definition()
+ #print sem_mesh.block_mat
+ #print sem_mesh.block_flag
+ #
+ sem_mesh.write(path=path_exporting_mesh_SPECFEM3D)
+ print 'END SPECFEM3D exporting process......'
+
if __name__ == '__main__':
- path='MESH/'
- export2SESAME(path)
-
-# call by:
-# import cubit2specfem3d
-# cubit2specfem3d.export2SESAME('MESH')
+ path='.'
+ export2SPECFEM3D(path)
Added: seismo/3D/SPECFEM3D/trunk/CUBIT/exportlib.py
===================================================================
--- seismo/3D/SPECFEM3D/trunk/CUBIT/exportlib.py (rev 0)
+++ seismo/3D/SPECFEM3D/trunk/CUBIT/exportlib.py 2013-09-17 16:14:58 UTC (rev 22796)
@@ -0,0 +1,667 @@
+#!/usr/bin/env python
+#############################################################################
+# exportlib.py
+# this file is part of GEOCUBIT #
+# #
+# Created by Emanuele Casarotti #
+# Copyright (c) 2008 Istituto Nazionale di Geofisica e Vulcanologia #
+# #
+#############################################################################
+# #
+# GEOCUBIT is free software: you can redistribute it and/or modify #
+# it under the terms of the GNU General Public License as published by #
+# the Free Software Foundation, either version 3 of the License, or #
+# (at your option) any later version. #
+# #
+# GEOCUBIT is distributed in the hope that it will be useful, #
+# but WITHOUT ANY WARRANTY; without even the implied warranty of #
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the #
+# GNU General Public License for more details. #
+# #
+# You should have received a copy of the GNU General Public License #
+# along with GEOCUBIT. If not, see <http://www.gnu.org/licenses/>. #
+# #
+#############################################################################
+#
+try:
+ import start as start
+ cubit = start.start_cubit()
+except:
+ try:
+ import cubit
+ except:
+ print 'error importing cubit, check if cubit is installed'
+ pass
+
+def refine_closecurve(block=1001,closed_filenames=None,acis=True):
+ from utilities import load_curves
+ from boundary_definition import build_block_side,define_surf
+ from mesh_volume import refine_inside_curve
+ #
+ #
+ curves=[]
+ if not isinstance(closed_filenames,list): closed_filenames=[closed_filenames]
+ for f in closed_filenames:
+ print f
+ if acis:
+ curves=curves+load_curves(f)
+ print curves
+ blist=list(cubit.get_block_id_list())
+ try:
+ blist.remove(1001)
+ except:
+ pass
+ try:
+ blist.remove(1002)
+ except:
+ pass
+ try:
+ blist.remove(1003)
+ except:
+ pass
+ try:
+ blist.remove(1004)
+ except:
+ pass
+ try:
+ blist.remove(1005)
+ except:
+ pass
+ try:
+ blist.remove(1006)
+ except:
+ pass
+ id_top=max(blist)
+ cmd='group "coi" add node in hex in block '+str(id_top)
+ cubit.cmd(cmd)
+ #
+ id_inside_arc=None
+ for c in map(int,curves[0].split()): #curves is a list of one string
+ c1001 = cubit.get_exodus_element_count(1001, "block")
+ c1002 = cubit.get_exodus_element_count(1002, "block")
+ c1003 = cubit.get_exodus_element_count(1003, "block")
+ c1004 = cubit.get_exodus_element_count(1004, "block")
+ c1005 = cubit.get_exodus_element_count(1005, "block")
+ c1006 = cubit.get_exodus_element_count(1006, "block")
+ #
+ refine_inside_curve(c,ntimes=1,depth=1,block=block,surface=False)
+ blist=list(cubit.get_block_id_list())
+ cmd='create mesh geometry hex all except hex in block all feature_angle 135'
+ cubit.cmd(cmd)
+ blist_after=list(cubit.get_block_id_list())
+ [blist_after.remove(x) for x in blist]
+ id_inside=max(blist_after)
+ cmd='group "coi" add node in hex in block '+str(id_inside)
+ cubit.cmd(cmd)
+ if id_inside_arc:
+ cmd='del block '+str(id_inside-1)
+ cubit.cmd(cmd)
+ cmd='block '+str(id_inside)+' name "refined"'
+ cubit.cmd(cmd)
+ id_inside_arc=id_inside
+ #
+ _,_,_,_,_,top_surf,bottom_surf,surf_xmin,surf_ymin,surf_xmax,surf_ymax=define_surf()
+ #
+ c1001_after = cubit.get_exodus_element_count(1001, "block")
+ c1002_after = cubit.get_exodus_element_count(1002, "block")
+ c1003_after = cubit.get_exodus_element_count(1003, "block")
+ c1004_after = cubit.get_exodus_element_count(1004, "block")
+ c1005_after = cubit.get_exodus_element_count(1005, "block")
+ c1006_after = cubit.get_exodus_element_count(1006, "block")
+ entity='face'
+ if c1001_after != c1001:
+ refname=entity+'_topo'
+ build_block_side(top_surf,refname,obj=entity,id_0=1001)
+ #
+ if c1002_after != c1002:
+ refname=entity+'_bottom'
+ build_block_side(bottom_surf,refname,obj=entity,id_0=1002)
+ #
+ if c1003_after != c1003:
+ refname=entity+'_abs_xmin'
+ build_block_side(surf_xmin,refname,obj=entity,id_0=1003)
+ #
+ if c1004_after != c1004:
+ refname=entity+'_abs_ymin'
+ build_block_side(surf_ymin,refname,obj=entity,id_0=1004)
+ #
+ if c1005_after != c1005:
+ refname=entity+'_abs_xmax'
+ build_block_side(surf_xmax,refname,obj=entity,id_0=1005)
+ #
+ if c1006_after != c1006:
+ refname=entity+'_abs_ymax'
+ build_block_side(surf_ymax,refname,obj=entity,id_0=1006)
+ #
+ cmd='disassociate mesh from volume all'
+ cubit.cmd(cmd)
+ cmd='group "coi" add node in face in block 1001 1002 1003 1004 1005 1006'
+ cubit.cmd(cmd)
+ cubit.cmd('del vol all')
+ cubit.cmd('group "removedouble" add hex all except hex in block all')
+ cubit.cmd('delete hex in removedouble')
+ cubit.cmd('delet group removedouble')
+ cmd='equivalence node in group coi tolerance 20'
+ cubit.cmd(cmd)
+ cmd='equivalence node all tolerance 10'
+ cubit.cmd(cmd)
+ cubit.cmd('del curve '+' '.join(str(x) for x in curves) )
+
+
+
+
+def collecting_merging(cpuxmin=0,cpuxmax=1,cpuymin=0,cpuymax=1,cpux=1,cpuy=1,cubfiles=False,ckbound_method1=False,ckbound_method2=False,merge_tolerance=None):
+ import glob
+ import re
+ #
+ rule_st=re.compile("(.+)_[0-9]+\.")
+ rule_ex=re.compile(".+_[0-9]+\.(.+)")
+ rule_int=re.compile(".+_([0-9]+)\.")
+ boundary_dict={}
+ ##
+ try:
+ from boundary_definition import check_bc, map_boundary
+ except:
+ pass
+ #
+ xmin,xmax,ymin,ymax,listfull=map_boundary(cpuxmin,cpuxmax,cpuymin,cpuymax,cpux,cpuy)
+ #
+ if cubfiles:
+ filenames=glob.glob(cubfiles)
+ try:
+ st = rule_st.findall(filenames[0])[0]
+ ex = rule_ex.findall(filenames[0])[0]
+ listflag=True
+ except:
+ ex=''
+ listflag=False
+
+ if ex == 'cub':
+ cubflag=True
+ else:
+ cubflag=False
+ list_int=[]
+ fs=[]
+ try:
+ for f in filenames:
+ i=int(rule_int.findall(f)[0])
+ list_int.append(i)
+ list_int.sort()
+ for i,ind in enumerate(list_int):
+ f=st+'_'+str(ind)+'.'+ex
+ fs.append(f)
+ except:
+ pass
+ if listflag:
+ filenames=fs
+ else:
+ pass
+ print len(filenames),filenames
+ else:
+ filenames=[]
+ ip=0
+ #
+ if len(filenames) > 0:
+ for filename in filenames[:]:
+ try:
+ ip=int(rule_int.findall(filename)[0])
+ if ip in listfull:
+ if cubflag:
+ cubit.cmd('import cubit "'+filename+'"')
+ else:
+ cubit.cmd('import mesh geometry "'+filename+'" block all use nodeset sideset feature_angle 135.00 linear merge')
+ boundary=check_bc(ip,xmin,xmax,ymin,ymax,cpux,cpuy,cpuxmin,cpuxmax,cpuymin,cpuymax)
+ boundary_dict[ip]=boundary
+ list_vol=list(cubit.parse_cubit_list('volume','all'))
+ for v in list_vol:
+ cubit.cmd("disassociate mesh from volume "+str(v))
+ command = "del vol "+str(v)
+ cubit.cmd(command)
+ except:
+ cubit.cmd('import mesh geometry "'+filename+'" block all use nodeset sideset feature_angle 135.00 linear merge')
+ ip=0
+ boundary=check_bc(ip,xmin,xmax,ymin,ymax,cpux,cpuy,cpuxmin,cpuxmax,cpuymin,cpuymax)
+ boundary_dict[ip]=boundary
+ list_vol=list(cubit.parse_cubit_list('volume','all'))
+ for v in list_vol:
+ cubit.cmd("disassociate mesh from volume "+str(v))
+ command = "del vol "+str(v)
+ cubit.cmd(command)
+ cubit.cmd('export mesh "tmp_collect_NOmerging.e" dimension 3 block all overwrite')
+ else:
+ boundary=check_bc(ip,xmin,xmax,ymin,ymax,cpux,cpuy,cpuxmin,cpuxmax,cpuymin,cpuymax)
+ #
+ #
+ block_list=cubit.get_block_id_list()
+ for block in block_list:
+ ty=cubit.get_block_element_type(block)
+ if ty == 'HEX8':
+ cubit.cmd('block '+str(block)+' name "vol'+str(block)+'"')
+ #
+ #
+ print 'chbound',ckbound_method1,ckbound_method2
+ if ckbound_method1 and not ckbound_method2:
+ if isinstance(merge_tolerance,list):
+ tol=merge_tolerance[0]
+ elif merge_tolerance:
+ tol=merge_tolerance
+ else:
+ tol=100000
+ #
+ idiag=None
+ #cubit.cmd('set info off')
+ #cubit.cmd('set journal off')
+ #cubit.cmd('set echo off')
+ ind=0
+ for ix in range(cpuxmin,cpuxmax):
+ for iy in range(cpuymin,cpuymax):
+ ind=ind+1
+ ip=iy*cpux+ix
+ print '******************* ',ip, ind,'/',len(listfull)
+ #
+ # ileft | ip
+ # --------------------
+ # idiag | idown
+ #
+ #
+ if ip not in xmin and ip not in ymin:
+ ileft=iy*cpux+ix-1
+ idown=(iy-1)*cpux+ix
+ idiag=idown-1
+ elif ip in xmin and ip in ymin:
+ ileft=ip
+ idown=ip
+ idiag=None
+ elif ip in xmin:
+ ileft=ip
+ idown=(iy-1)*cpux+ix
+ idiag=idown
+ elif ip in ymin:
+ ileft=iy*cpux+ix-1
+ idown=ip
+ idiag=ileft
+ #
+ if ip != idown:
+ nup=boundary_dict[ip]['nodes_surf_ymin']
+ ndow=boundary_dict[idown]['nodes_surf_ymax']
+ merge_node(nup,ndow)
+
+ if idiag != idown:
+ if ip in ymax and ip not in xmin:
+ nlu=boundary_dict[ip]['node_curve_xminymax'] #node in curve chunck left up... r u
+ nru=boundary_dict[ileft]['node_curve_xmaxymax']
+ merge_node(nlu,nru)
+ if ip in xmax:
+ nrd=boundary_dict[ip]['node_curve_xmaxymin'] #node in curve chunck left up... r u
+ nru=boundary_dict[idown]['node_curve_xmaxymax']
+ merge_node(nrd,nru)
+ nru=boundary_dict[ip]['node_curve_xminymin'] #node in curve chunck right up... r u
+ nrd=boundary_dict[idown]['node_curve_xminymax']
+ nld=boundary_dict[idiag]['node_curve_xmaxymax']
+ nlu=boundary_dict[ileft]['node_curve_xmaxymin']
+ merge_node_4(nru,nrd,nld,nlu)
+ elif ip in xmin:
+ nru=boundary_dict[ip]['node_curve_xminymin'] #node in curve chunck right up... r u
+ nrd=boundary_dict[idown]['node_curve_xminymax']
+ merge_node(nrd,nru)
+
+ #
+ if ip != ileft:
+ nright=boundary_dict[ip]['nodes_surf_xmin']
+ nleft=boundary_dict[ileft]['nodes_surf_xmax']
+ merge_node(nright,nleft)
+ #
+ #
+ if ip in ymin:
+ nrd=boundary_dict[ip]['node_curve_xminymin'] #node in curve chunck right down... r u
+ nld=boundary_dict[ileft]['node_curve_xmaxymin']
+ merge_node(nrd,nld)
+ if ip in ymax:
+ nru=boundary_dict[ip]['node_curve_xminymax'] #node in curve chunck right up... r u
+ nlu=boundary_dict[ileft]['node_curve_xmaxymax']
+ merge_node(nlu,nru)
+
+ cubit.cmd('set info on')
+ cubit.cmd('set echo on')
+ cubit.cmd('set journal on')
+
+
+ #
+ #
+ cmd='group "negativejac" add quality hex all Jacobian high'
+ cubit.cmd(cmd)
+ 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'
+ elif ckbound_method2 and not ckbound_method1:
+ if isinstance(merge_tolerance,list):
+ tol=merge_tolerance[0]
+ elif merge_tolerance:
+ tol=merge_tolerance
+ else:
+ tol=100000
+ #
+ idiag=None
+ for ix in range(cpuxmin,cpuxmax):
+ for iy in range(cpuymin,cpuymax):
+ ip=iy*cpux+ix
+ print '******************* ',ip
+ #
+ # ileft | ip
+ # --------------------
+ # idiag | idown
+ #
+ #
+ if ip not in xmin and ip not in ymin:
+ ileft=iy*cpux+ix-1
+ idown=(iy-1)*cpux+ix
+ idiag=idown-1
+ elif ip in xmin and ip in ymin:
+ ileft=ip
+ idown=ip
+ elif ip in xmin:
+ ileft=ip
+ idown=(iy-1)*cpux+ix
+ idiag=idown
+ elif ip in ymin:
+ ileft=iy*cpux+ix-1
+ idown=ip
+ idiag=ileft
+ #
+ #
+ if ip != idown:
+ nup=boundary_dict[ip]['nodes_surf_ymin']
+ ndow=boundary_dict[idown]['nodes_surf_ymax']
+ for n1,n2 in zip(nup,ndow):
+ cubit.cmd('equivalence node '+str(n1)+' '+str(n2)+' tolerance '+str(tol))
+ if idiag != idown:
+ if ip in ymax and ip not in xmin:
+ nlu=boundary_dict[ip]['node_curve_xminymax'] #node in curve chunck left up... r u
+ nru=boundary_dict[ileft]['node_curve_xmaxymax']
+ for n in zip(nlu,nru):
+ cubit.cmd('equivalence node '+' '.join(str(x) for x in n)+' tolerance '+str(tol))
+ nru=boundary_dict[ip]['node_curve_xminymin'] #node in curve chunck right up... r u
+ nrd=boundary_dict[idown]['node_curve_xminymax']
+ nld=boundary_dict[idiag]['node_curve_xmaxymax']
+ nlu=boundary_dict[ileft]['node_curve_xmaxymin']
+ for n in zip(nru,nrd,nlu,nld):
+ cubit.cmd('equivalence node '+' '.join(str(x) for x in n)+' tolerance '+str(tol))
+ elif ip in xmin:
+ nru=boundary_dict[ip]['node_curve_xminymin'] #node in curve chunck right up... r u
+ nrd=boundary_dict[idown]['node_curve_xminymax']
+ for n in zip(nru,nrd):
+ cubit.cmd('equivalence node '+' '.join(str(x) for x in n)+' tolerance '+str(tol))
+ #
+ #
+ if ip != ileft:
+ nright=boundary_dict[ip]['nodes_surf_xmin']
+ nleft=boundary_dict[ileft]['nodes_surf_xmax']
+ for n1,n2 in zip(nleft,nright):
+ cubit.cmd('equivalence node '+str(n1)+' '+str(n2)+' tolerance '+str(tol))
+ #
+ #
+ if ip in ymin:
+ nrd=boundary_dict[ip]['node_curve_xminymin'] #node in curve chunck right down... r u
+ nld=boundary_dict[ileft]['node_curve_xmaxymin']
+ for n in zip(nrd,nld):
+ cubit.cmd('equivalence node '+' '.join(str(x) for x in n)+' tolerance '+str(tol))
+ if ip in ymax:
+ nru=boundary_dict[ip]['node_curve_xminymax'] #node in curve chunck right up... r u
+ nlu=boundary_dict[ileft]['node_curve_xmaxymax']
+ for n in zip(nru,nlu):
+ cubit.cmd('equivalence node '+' '.join(str(x) for x in n)+' tolerance '+str(tol))
+ #
+ #
+ cmd='topology check coincident node face all tolerance '+str(tol*2)+' nodraw brief result group "checkcoinc"'
+ cubit.silent_cmd(cmd)
+ group_id_1=cubit.get_id_from_name("checkcoinc")
+ 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'
+ import sys
+ sys.exit()
+ cmd='group "negativejac" add quality hex all Jacobian high'
+ cubit.cmd(cmd)
+ 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'
+ elif ckbound_method1 and ckbound_method2:
+ block_list=cubit.get_block_id_list()
+ i=-1
+ for block in block_list:
+ ty=cubit.get_block_element_type(block)
+ if ty == 'HEX8':
+ i=i+1
+ if isinstance(merge_tolerance,list):
+ try:
+ tol=merge_tolerance[i]
+ except:
+ tol=merge_tolerance[-1]
+ elif merge_tolerance:
+ tol=merge_tolerance
+ else:
+ tol=1
+ cmd='topology check coincident node face in hex in block '+str(block)+' tolerance '+str(tol)+' nodraw brief result group "b'+str(block)+'"'
+ cubit.cmd(cmd)
+ print cmd
+ cmd='equivalence node in group b'+str(block)+' tolerance '+str(tol)
+ cubit.cmd(cmd)
+ print cmd
+ if isinstance(merge_tolerance,list):
+ tol=max(merge_tolerance)
+ elif merge_tolerance:
+ tol=merge_tolerance
+ else:
+ tol=1
+ #
+ #
+ cmd='topology check coincident node face all tolerance '+str(tol)+' nodraw brief result group "checkcoinc"'
+ cubit.silent_cmd(cmd)
+ group_id_1=cubit.get_id_from_name("checkcoinc")
+ 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'
+ import sys
+ sys.exit()
+ cmd='group "negativejac" add quality hex all Jacobian high'
+ cubit.silent_cmd(cmd)
+ 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'
+
+def collect(cpuxmin=0,cpuxmax=1,cpuymin=0,cpuymax=1,cpux=1,cpuy=1,cubfiles=False,ckbound_method1=False,ckbound_method2=False,merge_tolerance=None,curverefining=False,outfilename='totalmesh_merged',qlog=False,export2SPECFEM3D=False,listblock=None,listflag=None,outdir='.'):
+ #
+ collecting_merging(cpuxmin,cpuxmax,cpuymin,cpuymax,cpux,cpuy,cubfiles=cubfiles,ckbound_method1=ckbound_method1,ckbound_method2=ckbound_method2,merge_tolerance=merge_tolerance)
+ #
+ if curverefining:
+ block=1001 #topography
+ refine_closecurve(block,curverefining,acis=True)
+ #
+ #
+ cubit.cmd('compress all')
+ command="export mesh '"+outfilename+".e' block all overwrite xml '"+outfilename+".xml'"
+ cubit.cmd(command)
+ outdir2='/'.join(x for x in outfilename.split('/')[:-1])
+ if outdir2 == '': outdir2='.'
+ f=open(outdir2+'/'+'blocks.dat','w')
+ blocks=cubit.get_block_id_list()
+ #
+ for block in blocks:
+ name=cubit.get_exodus_entity_name('block',block)
+ element_count = cubit.get_exodus_element_count(block, "block")
+ nattrib=cubit.get_block_attribute_count(block)
+ attr=[cubit.get_block_attribute_value(block,x) for x in range(0,nattrib)]
+ ty=cubit.get_block_element_type(block)
+ f.write(str(block)+' ; '+name+' ; nattr '+str(nattrib)+' ; '+' '.join(str(x) for x in attr)+' ; '+ty+' '+str(element_count)+'\n')
+ f.close()
+ #
+ #
+ cubit.cmd('set info echo journ off')
+ cmd='del group all'
+ cubit.silent_cmd(cmd)
+ cubit.cmd('set info echo journ on')
+ #
+ command = "save as '"+outfilename+".cub' overwrite"
+ cubit.cmd(command)
+ #
+ print 'end meshing'
+ #
+ #
+ if qlog:
+ print '\n\nQUALITY CHECK.... ***************\n\n'
+ import quality_log
+ tq=open(outfilename+'.quality','w')
+ max_skewness,min_length=quality_log.quality_log(tq)
+ #
+ #
+ #
+ if export2SPECFEM3D:
+ e2SEM(files=False,listblock=listblock,listflag=listflag,outdir=outdir)
+
+def e2SEM(files=False,listblock=None,listflag=None,outdir='.'):
+ import glob
+ if files:
+ filenames=glob.glob(files)
+ for f in filenames:
+ 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
+ if listblock and listflag:
+ pass
+ else:
+ listblock=[]
+ listflag=[]
+ block_list=list(cubit.get_block_id_list())
+ for block in block_list:
+ ty=cubit.get_block_element_type(block)
+ if 'HEX' in ty:
+ listblock.append(block)
+ #listflag.append(block)
+ listflag=range(1,len(listflag)+1)
+ #
+ for ib,iflag in zip(listblock,listflag):
+ cubit.cmd("block "+str(ib)+" attribute count 1")
+ cubit.cmd("block "+str(ib)+" attribute index 1 "+ str(iflag) )
+ #
+ import cubit2specfem3d
+ cubit2specfem3d.export2SPECFEM3D(outdir)
+
+def invert_dict(d):
+ inv = {}
+ for k,v in d.iteritems():
+ keys = inv.setdefault(v, [])
+ keys.append(k)
+ return inv
+
+
+
+def prepare_equivalence(nodes1,nodes2):
+ cubit.cmd('set info off')
+ cubit.cmd('set echo off')
+ cubit.cmd('set journal off')
+ length={}
+ for ns in zip(nodes1,nodes2):
+ cmd='group "tmpn" add edge in node '+' '.join(str(n) for n in ns )
+ cubit.cmd(cmd)
+ ge=cubit.get_id_from_name("tmpn")
+ e1=cubit.get_group_edges(ge)
+ lengthmin=1e9
+ for e in e1:
+ lengthmin=min(lengthmin,cubit.get_mesh_edge_length(e))
+ length[ns]=lengthmin*.5
+ cubit.cmd('delete group '+str(ge))
+ minvalue=min(length.values())
+ maxvalue=max(length.values())
+ print 'min lentgh: ',minvalue,'max lentgh: ',maxvalue
+ nbin= int((maxvalue/minvalue)/2.)+1
+ factor=(maxvalue-minvalue)/nbin
+ dic_new={}
+ for k in length.keys():
+ dic_new[k]=int((length[k]-minvalue)/factor)
+ inv_length=invert_dict(dic_new)
+ print inv_length.keys(),factor,minvalue
+ ks=inv_length.keys()
+ ks.sort()
+ for k in range(0,len(inv_length.keys())-1):
+ inv_length[ks[k]]=inv_length[ks[k]]+inv_length[ks[k+1]]
+ cubit.cmd('set info on')
+ cubit.cmd('set echo on')
+ cubit.cmd('set journal on')
+ return factor,minvalue,inv_length
+
+
+def merge_node(n1,n2):
+ factor,minvalue,inv_length=prepare_equivalence(n1,n2)
+
+ for k in inv_length.keys()[:-1]:
+ if len(inv_length[k]) > 0:
+ cmd='equivalence node '+' '.join(' '.join(str(n) for n in x) for x in inv_length[k])+' tolerance '+str(k*factor+minvalue/2.)
+ cubit.cmd(cmd)
+ print 'equivalence '+str(len(inv_length[k]))+' couples of nodes - tolerance '+str(k*factor+minvalue/2.)
+
+
+def prepare_equivalence_4(nodes1,nodes2,nodes3,nodes4):
+ cubit.cmd('set info off')
+ cubit.cmd('set echo off')
+ cubit.cmd('set journal off')
+ length={}
+ for ns in zip(nodes1,nodes2,nodes3,nodes4):
+ cmd='group "tmpn" add edge in node '+' '.join(str(n) for n in ns )
+ cubit.cmd(cmd)
+ ge=cubit.get_id_from_name("tmpn")
+ e1=cubit.get_group_edges(ge)
+ lengthmin=1e9
+ for e in e1:
+ lengthmin=min(lengthmin,cubit.get_mesh_edge_length(e))
+ length[ns]=lengthmin*.5
+ cubit.cmd('delete group '+str(ge))
+ try:
+ minvalue=min(length.values())
+ maxvalue=max(length.values())
+ except:
+ print nodes1,nodes2,nodes3,nodes4
+ print 'edges ', e1
+ minvalue=100.
+ maxvalue=2000.
+ print 'min lentgh: ',minvalue,'max lentgh: ',maxvalue
+ nbin= int((maxvalue/minvalue)/2.)+1
+ factor=(maxvalue-minvalue)/nbin
+ dic_new={}
+ for k in length.keys():
+ dic_new[k]=int((length[k]-minvalue)/factor)
+ inv_length=invert_dict(dic_new)
+ print inv_length.keys(),factor,minvalue
+ ks=inv_length.keys()
+ ks.sort()
+ for k in range(0,len(inv_length.keys())-1):
+ inv_length[ks[k]]=inv_length[ks[k]]+inv_length[ks[k+1]]
+ cubit.cmd('set info on')
+ cubit.cmd('set echo on')
+ cubit.cmd('set journal on')
+ return factor,minvalue,inv_length
+
+
+def merge_node_4(n1,n2,n3,n4):
+ factor,minvalue,inv_length=prepare_equivalence_4(n1,n2,n3,n4)
+
+ for k in inv_length.keys()[:-1]:
+ if len(inv_length[k]) > 0:
+ cmd='equivalence node '+' '.join(' '.join(str(n) for n in x) for x in inv_length[k])+' tolerance '+str(k*factor+minvalue/2.)
+ cubit.cmd(cmd)
+ print 'equivalence '+str(len(inv_length[k]))+' couples of nodes - tolerance '+str(k*factor+minvalue/2.)
+
+
+
+
+
+
\ No newline at end of file
Added: seismo/3D/SPECFEM3D/trunk/CUBIT/hex_metric.py
===================================================================
--- seismo/3D/SPECFEM3D/trunk/CUBIT/hex_metric.py (rev 0)
+++ seismo/3D/SPECFEM3D/trunk/CUBIT/hex_metric.py 2013-09-17 16:14:58 UTC (rev 22796)
@@ -0,0 +1,772 @@
+#############################################################################
+# hex_metric.py
+# this file is part of GEOCUBIT #
+# #
+# Created by Emanuele Casarotti #
+# Copyright (c) 2008 Istituto Nazionale di Geofisica e Vulcanologia #
+# #
+#############################################################################
+# #
+# GEOCUBIT is free software: you can redistribute it and/or modify #
+# it under the terms of the GNU General Public License as published by #
+# the Free Software Foundation, either version 3 of the License, or #
+# (at your option) any later version. #
+# #
+# GEOCUBIT is distributed in the hope that it will be useful, #
+# but WITHOUT ANY WARRANTY; without even the implied warranty of #
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the #
+# GNU General Public License for more details. #
+# #
+# You should have received a copy of the GNU General Public License #
+# along with GEOCUBIT. If not, see <http://www.gnu.org/licenses/>. #
+# #
+#############################################################################
+##mesh=SEM_metric_3D
+##mesh.check_metric()
+##print mesh
+#
+#
+import math
+
+try:
+ import start as start
+ cubit = start.start_cubit()
+except:
+ try:
+ import cubit
+ except:
+ print 'error importing cubit, check if cubit is installed'
+ pass
+
+class SEM_metric_3D(object):
+ def __init__(self,list_hex=None,volume=None):
+ super(SEM_metric_3D, self).__init__()
+ self.list_hex=list_hex
+ self.volume=volume
+ self.skew_hystogram=None
+ self.max_skewness=None
+ self.max_angle=None
+ self.min_angle=None
+ self.min_edge_length=None
+ self.cubit_skew=0
+ self.max_valence=None
+ self.node_with_max_valence=None
+ self.valence_threshold=7
+ self.resolution=.05
+ self.resolution_edge_length=.1
+ self.nbin=10
+ self.spatial_unit=1
+ self.gll=0.17267316464601141
+ self.dt=None
+ 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
+ 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,')'
+ 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
+ tot=0
+ for i in self.skew_hystogram.values():
+ tot=tot+len(i)
+ #k=self.skew_hystogram.keys()
+ #k.sort()
+ factor=self.max_skewness/self.nbin
+ for i in range(0,self.nbin+1):
+ if self.skew_hystogram.has_key(i):
+ if (i+1)*factor <= 1:
+ print i,' [',i*factor,'->',(i+1)*factor,'[ : ',len(self.skew_hystogram[i]),'/',tot,' hexes (',len(self.skew_hystogram[i])/float(tot)*100.,'%)'
+ else:
+ if (i+1)*factor <= 1:
+ print i,' [',i*factor,'->',(i+1)*factor,'[ : ',0,'/',tot,' hexes (0%)'
+ print
+ ###############################################
+ if self.min_edge_length is not None:
+ self.edgemin_hystogram=self.hyst(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
+ if len(self.hex_min_edge_length) <= 30:
+ 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.'
+ if len(self.hex_max_edge_length) <= 30:
+ 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
+ factor=(self.max_edge_length-self.min_edge_length)/self.nbin
+ print 'minimum edge length'
+ tot=0
+ for i in self.edgemin_hystogram.values():
+ tot=tot+len(i)
+ #k=self.edgemin_hystogram.keys()
+ #k.sort()
+ for i in range(0,self.nbin+1):
+ if self.edgemin_hystogram.has_key(i):
+ print i,' [',i*factor+self.min_edge_length,'->',(i+1)*factor+self.min_edge_length,'[ : ',len(self.edgemin_hystogram[i]),'/',tot,' hexes (',len(self.edgemin_hystogram[i])/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'
+ tot=0
+ for i in self.edgemax_hystogram.values():
+ tot=tot+len(i)
+ #k=self.edgemax_hystogram.keys()
+ #k.sort()
+ for i in range(0,self.nbin+1):
+ if self.edgemax_hystogram.has_key(i):
+ print i,' [',i*factor+self.min_edge_length,'->',(i+1)*factor+self.min_edge_length,'[ : ',len(self.edgemax_hystogram[i]),'/',tot,' hexes (',len(self.edgemax_hystogram[i])/float(tot)*100.,'%)'
+ else:
+ 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
+ try:
+ return str(len(self.list_hex))+' hexes checked'
+ except:
+ return 'please performe a metric check: ex ... object.check_metric()'
+ #
+ #
+ def invert_dict(self,d):
+ inv = {}
+ for k,v in d.iteritems():
+ keys = inv.setdefault(v, [])
+ keys.append(k)
+ return inv
+ def hex_metric(self,h):
+ nodes=cubit.get_connectivity('Hex',h)
+ equiangle_skewness=None
+ min_angle=float('infinity')
+ edge_length_min=float('infinity')
+ max_angle=None
+ edge_length_max=None
+ #loopnode=[(i,j) for i in range(8) for j in range(i+1,8)]
+ if len(nodes) == 4:
+ faces=[[0,1,2,3]]
+ elif len(nodes) == 8:
+ 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'
+ return None,None,None
+ x=[]
+ y=[]
+ z=[]
+ for i in nodes:
+ n=cubit.get_nodal_coordinates(i)
+ x.append(n[0])
+ y.append(n[1])
+ z.append(n[2])
+ for face in faces:
+ for i in range(-1,3):
+ vx1=x[face[i-1]]-x[face[i]]
+ vy1=y[face[i-1]]-y[face[i]]
+ vz1=z[face[i-1]]-z[face[i]]
+ #
+ vx2=x[face[i+1]]-x[face[i]]
+ vy2=y[face[i+1]]-y[face[i]]
+ vz2=z[face[i+1]]-z[face[i]]
+ #
+ norm1=math.sqrt(vx1*vx1+vy1*vy1+vz1*vz1)
+ norm2=math.sqrt(vx2*vx2+vy2*vy2+vz2*vz2)
+ #
+ if norm1 == 0 or norm2 == 0:
+ print 'degenerated mesh, 0 length edge'
+ import sys
+ sys.exit()
+ angle=math.acos((vx1*vx2+vy1*vy2+vz1*vz2)/(norm1*norm2))
+ #
+ equiangle_skewness=max(equiangle_skewness,math.fabs(2.*angle-math.pi)/math.pi)
+ min_angle=min(min_angle,angle*180./math.pi)
+ max_angle=max(max_angle,angle*180./math.pi)
+ edge_length_min=min(edge_length_min,norm1)
+ edge_length_max=max(edge_length_max,norm1)
+ return round(equiangle_skewness,5),min_angle,max_angle,round(edge_length_min,5),round(edge_length_max,5)
+ def show(self,minvalue,maxvalue,dic):
+ hexlist=[]
+ for k in dic.keys():
+ if minvalue <= dic[k] <= maxvalue: hexlist.extend([k])
+ #command = "draw hex "+str(hexlist)
+ #command = command.replace("["," ").replace("]"," ").replace("("," ").replace(")"," ")
+ #cubit.cmd(command)
+ return hexlist
+ def hyst(self,minvalue,maxvalue,dic):
+ if maxvalue != minvalue:
+ factor=(maxvalue-minvalue)/self.nbin
+ else:
+ factor=1
+ dic_new={}
+ for k in dic.keys():
+ dic_new[k]=int((dic[k]-minvalue)/factor)
+ inv_dic=self.invert_dict(dic_new)
+ return inv_dic
+ #
+ #
+ #
+ #
+ #
+ ############################################################################################
+ def pick_hex(self,list_hex=None,volume=None):
+ if self.cubit_skew > 0:
+ command = "del group skew_top"
+ cubit.cmd(command)
+ command = "group 'skew_top' add quality volume all skew low "+str(self.cubit_skew)
+ cubit.silent_cmd(command)
+ group=cubit.get_id_from_name("skew_top")
+ self.list_hex=cubit.get_group_hexes(group)
+ elif list_hex is not None:
+ self.list_hex=list_hex
+ elif volume is not None:
+ command = "group 'hextmp' add hex in volume "+str(volume)
+ cubit.silent_cmd(command)
+ group=cubit.get_id_from_name("hextmp")
+ self.list_hex=cubit.get_group_hexes(group)
+ command = "del group hextmp"
+ cubit.silent_cmd(command)
+ elif self.volume is not None:
+ command = "group 'hextmp' add hex in volume "+str(self.volume)
+ cubit.silent_cmd(command)
+ group=cubit.get_id_from_name("hextmp")
+ self.list_hex=cubit.get_group_hexes(group)
+ command = "del group hextmp"
+ 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'
+ #
+ #
+ ############################################################################################
+ #SKEWNESS
+ def check_skew(self,list_hex=None,volume=None):
+ self.pick_hex(list_hex=list_hex,volume=volume)
+ global_max_skew=None
+ global_max_angle=None
+ global_min_angle=float('infinity')
+ tmp=int(1./self.resolution)+1
+ skew_hyst={}
+ for i in range(0,tmp):
+ if i*self.resolution < 1:
+ skew_hyst[i]=[]
+ h_global_max_skew=[]
+ h_global_max_angle=[]
+ h_global_min_angle=[]
+ for h in self.list_hex:
+ equiangle_skewness,min_angle,max_angle,tmp,tmp2=self.hex_metric(h)
+ skew_hyst[h]=equiangle_skewness
+ #
+ if equiangle_skewness > global_max_skew:
+ global_max_skew=equiangle_skewness
+ h_global_max_shew=[]
+ h_global_max_shew.append(h)
+ global_min_angle=min_angle
+ global_max_angle=max_angle
+ elif equiangle_skewness == global_max_skew:
+ h_global_max_shew.append(h)
+ s=h_global_max_shew
+ self.skew_hyst=skew_hyst
+ self.max_skewness=global_max_skew
+ self.max_angle=global_max_angle
+ self.min_angle=global_min_angle
+ self.hex_max_skewness=s
+ def show_skew(self,low=None,high=None):
+ if low is None and high is not None:
+ hexlist=self.show(low,self.max_skewness,self.skew_hyst)
+ elif low is not None and high is None:
+ hexlist=self.show(0,high,self.skew_hyst)
+ if low is not None and high is not None:
+ hexlist=self.show(low,high,self.skew_hyst)
+ else:
+ hexlist=self.hex_min_edge_length
+ command = "draw hex "+str(hexlist)
+ command = command.replace("["," ").replace("]"," ").replace("("," ").replace(")"," ")
+ cubit.silent_cmd(command)
+ def list_skew(self,low=None,high=None):
+ if low is None and high is not None:
+ hexlist=self.show(low,self.max_skewness,self.skew_hyst)
+ elif low is not None and high is None:
+ hexlist=self.show(0,high,self.skew_hyst)
+ if low is not None and high is not None:
+ hexlist=self.show(low,high,self.skew_hyst)
+ else:
+ hexlist=self.hex_min_edge_length
+ return hexlist
+ def highlight_skew(self,low=None,high=None):
+ if low is None and high is not None:
+ hexlist=self.show(low,self.max_skewness,self.skew_hyst)
+ elif low is not None and high is None:
+ hexlist=self.show(0,high,self.skew_hyst)
+ if low is not None and high is not None:
+ hexlist=self.show(low,high,self.skew_hyst)
+ else:
+ hexlist=self.hex_min_edge_length
+ command = "highlight hex "+str(hexlist)
+ command = command.replace("["," ").replace("]"," ").replace("("," ").replace(")"," ")
+ cubit.silent_cmd(command)
+ #
+ #
+ #
+ #############################################################################################
+ #HEX VALENCE
+ def check_valence(self):
+ #usage: hex_valence()
+ #
+ list_hex=cubit.parse_cubit_list('hex','all')
+ list_node=cubit.parse_cubit_list('node','all')
+ lookup=dict(zip(list_node,[0]*len(list_node)))
+ for h in list_hex:
+ n=cubit.get_connectivity('Hex',h)
+ for i in n:
+ if lookup.has_key(i): lookup[i]=lookup[i]+1
+ inv_lookup=self.invert_dict(lookup)
+ #
+ vmax=max(inv_lookup.keys())
+ nmax=inv_lookup[max(inv_lookup.keys())]
+ 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])))
+ self.valence_summary=inv_lookup
+ self.max_valence=vmax
+ self.node_with_max_valence=nmax
+ #
+ #
+ ###########################################################################################
+ #EDGE LENGTH
+ def check_edge_length(self,list_hex=None,volume=None):
+ self.pick_hex(list_hex=list_hex,volume=volume)
+ global_min_edge_length=float('infinity')
+ h_global_min_edge_length=[]
+ global_max_edge_length=None
+ h_global_max_edge_length=[]
+ edgemin_hyst={}
+ edgemax_hyst={}
+ for h in self.list_hex:
+ tmp,tmp2,tmp3,edge_length_min,edge_length_max=self.hex_metric(h)
+ edgemax_hyst[h]=edge_length_max
+ edgemin_hyst[h]=edge_length_min
+ if edge_length_min < global_min_edge_length:
+ global_min_edge_length=edge_length_min
+ h_global_min_edge_length=[]
+ h_global_min_edge_length.append(h)
+ elif edge_length_min == global_min_edge_length:
+ h_global_min_edge_length.append(h)
+ if edge_length_max > global_max_edge_length:
+ global_max_edge_length=edge_length_max
+ h_global_max_edge_length=[]
+ h_global_max_edge_length.append(h)
+ elif edge_length_max == global_max_edge_length:
+ h_global_max_edge_length.append(h)
+ self.min_edge_length=global_min_edge_length
+ self.hex_min_edge_length=h_global_min_edge_length
+ self.max_edge_length=global_max_edge_length
+ self.hex_max_edge_length=h_global_max_edge_length
+ self.edgemin_hyst=edgemin_hyst
+ self.edgemax_hyst=edgemax_hyst
+ def show_edgemin(self,low=None,high=None):
+ if low is None and high is not None:
+ hexlist=self.show(low,self.max_edge_length,self.edgemin_hyst)
+ elif low is not None and high is None:
+ hexlist=self.show(self.min_edge_length,high,self.edgemin_hyst)
+ if low is not None and high is not None:
+ hexlist=self.show(low,high,self.edgemin_hyst)
+ else:
+ hexlist=self.hex_min_edge_length
+ command = "draw hex "+str(hexlist)
+ command = command.replace("["," ").replace("]"," ").replace("("," ").replace(")"," ")
+ cubit.cmd(command)
+ def show_edgemax(self,low=None,high=None):
+ if low is None and high is not None:
+ hexlist=self.show(low,self.max_edge_length,self.edgemax_hyst)
+ elif low is not None and high is None:
+ hexlist=self.show(self.min_edge_length,high,self.edgemax_hyst)
+ if low is not None and high is not None:
+ hexlist=self.show(low,high,self.edgemax_hyst)
+ else:
+ hexlist=self.hex_max_edge_length
+ command = "draw hex "+str(hexlist)
+ command = command.replace("["," ").replace("]"," ").replace("("," ").replace(")"," ")
+ cubit.cmd(command)
+ def list_edgemax(self,low=None,high=None):
+ if low is None and high is not None:
+ hexlist=self.show(low,self.max_edge_length,self.edgemax_hyst)
+ elif low is not None and high is None:
+ hexlist=self.show(self.min_edge_length,high,self.edgemax_hyst)
+ if low is not None and high is not None:
+ hexlist=self.show(low,high,self.edgemax_hyst)
+ else:
+ hexlist=self.hex_max_edge_length
+ return hexlist
+ def list_edgemin(self,low=None,high=None):
+ if low is None and high is not None:
+ hexlist=self.show(low,self.max_edge_length,self.edgemin_hyst)
+ elif low is not None and high is None:
+ hexlist=self.show(self.min_edge_length,high,self.edgemin_hyst)
+ if low is not None and high is not None:
+ hexlist=self.show(low,high,self.edgemin_hyst)
+ else:
+ hexlist=self.hex_min_edge_length
+ return hexlist
+ def highlight_edgemin(self,low=None,high=None):
+ if low is None and high is not None:
+ hexlist=self.show(low,self.max_edge_length,self.edgemin_hyst)
+ elif low is not None and high is None:
+ hexlist=self.show(self.min_edge_length,high,self.edgemin_hyst)
+ if low is not None and high is not None:
+ hexlist=self.show(low,high,self.edgemin_hyst)
+ else:
+ hexlist=self.hex_min_edge_length
+ command = "highlight hex "+str(hexlist)
+ command = command.replace("["," ").replace("]"," ").replace("("," ").replace(")"," ")
+ cubit.cmd(command)
+ def highlight_edgemax(self,low=None,high=None):
+ if low is None and high is not None:
+ hexlist=self.show(low,self.max_edge_length,self.edgemax_hyst)
+ elif low is not None and high is None:
+ hexlist=self.show(self.min_edge_length,high,self.edgemax_hyst)
+ if low is not None and high is not None:
+ hexlist=self.show(low,high,self.edgemax_hyst)
+ else:
+ hexlist=self.hex_max_edge_length
+ command = "highlight hex "+str(hexlist)
+ command = command.replace("["," ").replace("]"," ").replace("("," ").replace(")"," ")
+ cubit.cmd(command) #
+ #
+ #
+ #
+ #
+ ############################################################################################
+ def check_metric(self,list_hex=None,volume=None):
+ self.pick_hex(list_hex=list_hex,volume=volume)
+ skew_hyst={}
+ edata={}
+ edgemin_hyst={}
+ edgemax_hyst={}
+ edgemaxdata={}
+ edgemindata={}
+ #
+ global_min_edge_length=float('infinity')
+ h_global_min_edge_length=[]
+ global_max_edge_length=None
+ h_global_max_edge_length=[]
+ #
+ global_max_skew=None
+ global_max_angle=None
+ global_min_angle=float('infinity')
+ h_global_max_skew=[]
+ h_global_max_angle=[]
+ h_global_min_angle=[]
+ #
+ s=0.
+ for h in self.list_hex:
+ equiangle_skewness,min_angle,max_angle,edge_length_min,edge_length_max=self.hex_metric(h)
+ skew_hyst[h]=equiangle_skewness
+ edgemax_hyst[h]=edge_length_max
+ edgemin_hyst[h]=edge_length_min
+ if edge_length_min < global_min_edge_length:
+ global_min_edge_length=edge_length_min
+ h_global_min_edge_length=[]
+ h_global_min_edge_length.append(h)
+ elif edge_length_min == global_min_edge_length:
+ h_global_min_edge_length.append(h)
+ if edge_length_max > global_max_edge_length:
+ global_max_edge_length=edge_length_max
+ h_global_max_edge_length=[]
+ h_global_max_edge_length.append(h)
+ elif edge_length_max == global_max_edge_length:
+ h_global_max_edge_length.append(h)
+ if equiangle_skewness > global_max_skew:
+ global_max_skew=equiangle_skewness
+ h_global_max_shew=[]
+ h_global_max_shew.append(h)
+ global_min_angle=min_angle
+ global_max_angle=max_angle
+ elif equiangle_skewness == global_max_skew:
+ h_global_max_shew.append(h)
+ s=h_global_max_shew
+ #
+ self.max_skewness=global_max_skew
+ self.max_angle=global_max_angle
+ self.min_angle=global_min_angle
+ self.hex_max_skewness=s
+ self.skew_hyst=skew_hyst
+ #
+ self.min_edge_length=global_min_edge_length
+ self.hex_min_edge_length=h_global_min_edge_length
+ self.max_edge_length=global_max_edge_length
+ self.hex_max_edge_length=h_global_max_edge_length
+ self.edgemin_hyst=edgemin_hyst
+ self.edgemax_hyst=edgemax_hyst
+ #
+
+class SEM_stability_3D(SEM_metric_3D):
+ def __init__(self,list_hex=None,volume=None):
+ super(SEM_metric_3D, self).__init__()
+ self.list_hex=list_hex
+ self.volume=volume
+ self.Ngll_per_wavelength=5
+ self.gllcoeff=.17
+ self.maxgllcoeff=0.5-.17
+ self.Cmax=.3
+ self.nbin=10
+ self.period_hyst=None
+ self.dt=None
+ self.cubit_skew =0
+ #
+ #
+ #
+ def __repr__(self):
+ print 'check mesh stability'
+ #
+ def check_simulation_parameter(self,list_hex=None,volume=None,tomofile=None,vp_static=None,vs_static=None):
+ self.pick_hex(list_hex=list_hex,volume=volume)
+ timestep_hyst={}
+ stability_hyst={}
+ global_min_timestep=float(1)
+ global_max_periodresolved=None
+ if tomofile:
+ self.read_tomo(tomofile)
+ #
+ #
+ for ind,h in enumerate(self.list_hex):
+ if ind%10000==0: 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
+ stability_hyst[h]=pmax_tmp
+ if dt_tmp < global_min_timestep:
+ global_min_timestep=dt_tmp
+ if pmax_tmp > global_max_periodresolved:
+ global_max_periodresolved=pmax_tmp
+ #
+ self.period=global_max_periodresolved
+ self.dt=global_min_timestep
+ self.dt_hyst=timestep_hyst
+ self.period_hyst=stability_hyst
+ #
+ def read_tomo(self,tomofile=None):
+ if tomofile:
+ print 'reading tomography file ',tomofile
+ import numpy
+ #xtomo,ytomo,ztomo,vp,vs,rho=numpy.loadtxt(tomofile,skiprows=4)
+ 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())
+ spacing_x, spacing_y, spacing_z =map(float,tf.readline().split())
+ nx, ny, nz =map(int,tf.readline().split())
+ vp_min, vp_max, vs_min, vs_max, rho_min, rho_max =map(float,tf.readline().split())
+ #
+ ind=0
+ import sys
+ xtomo,ytomo,ztomo,vp,vs=[],[],[],[],[]
+ while ind<nx*ny*nz:
+ ind=ind+1
+ if ind%100000==0: sys.stdout.write("reading progress: %i/%i \r" % (ind,(nx*ny*nz)) )
+ x,y,z,v1,v2,_=map(float,tf.readline().split())
+ xtomo.append(x)
+ ytomo.append(y)
+ ztomo.append(z)
+ vp.append(v1)
+ vs.append(v2)
+ tf.close()
+ 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
+ self.spacing_x, self.spacing_y, self.spacing_z = spacing_x, spacing_y, spacing_z
+ self.nx, self.ny, self.nz = nx, ny, nz
+ self.vp_min, self.vp_max, self.vs_min, self.vs_max = vp_min, vp_max, vs_min, vs_max
+ self.xtomo,self.ytomo,self.ztomo,self.vp,self.vs=xtomo,ytomo,ztomo,vp,vs
+ else:
+ print 'no tomofile!!!!'
+ #
+ #
+ #
+ def tomo(self,x,y,z,vp_static=None,vs_static=None):
+ if not vp_static:
+ spac_x = (x - self.orig_x) / self.spacing_x
+ spac_y = (y - self.orig_y) / self.spacing_y
+ spac_z = (z - self.orig_z) / self.spacing_z
+ #
+ ix = int(spac_x)
+ iy = int(spac_y)
+ iz = int(spac_z)
+ #
+ gamma_interp_x = spac_x - float(ix)
+ gamma_interp_y = spac_y - float(iy)
+ #
+ NX=self.nx
+ NY=self.ny
+ NZ=self.nz
+ if(ix < 0):
+ ix = 0
+ gamma_interp_x = 0.
+ if(ix > NX-2):
+ ix = NX-2
+ gamma_interp_x = 1.
+ if(iy < 0):
+ iy = 0
+ gamma_interp_y = 0.
+ if(iy > NY-2):
+ iy = NY-2
+ gamma_interp_y = 1.
+ if(iz < 0):
+ iz = 0
+ if(iz > NZ-2):
+ iz = NZ-2
+ #
+ p0 = ix+iy*NX+iz*(NX*NY)
+ p1 = (ix+1)+iy*NX+iz*(NX*NY)
+ p2 = (ix+1)+(iy+1)*NX+iz*(NX*NY)
+ p3 = ix+(iy+1)*NX+iz*(NX*NY)
+ p4 = ix+iy*NX+(iz+1)*(NX*NY)
+ p5 = (ix+1)+iy*NX+(iz+1)*(NX*NY)
+ p6 = (ix+1)+(iy+1)*NX+(iz+1)*(NX*NY)
+ p7 = ix+(iy+1)*NX+(iz+1)*(NX*NY)
+ #
+ if self.ztomo[p4]==self.ztomo[p0]:
+ gamma_interp_z1 = 1
+ else:
+ gamma_interp_z1 = (z-self.ztomo[p0])/(self.ztomo[p4]-self.ztomo[p0])
+ if(gamma_interp_z1 > 1.): gamma_interp_z1 = 1.
+ if(gamma_interp_z1 < 0.): gamma_interp_z1 = 0.
+ #
+ if self.ztomo[p5]==self.ztomo[p1]:
+ gamma_interp_z2 = 1
+ else:
+ gamma_interp_z2 = (z-self.ztomo[p1])/(self.ztomo[p5]-self.ztomo[p1])
+ if(gamma_interp_z2 > 1.): gamma_interp_z2 = 1.
+ if(gamma_interp_z2 < 0.): gamma_interp_z2 = 0.
+ #
+ if self.ztomo[p6]==self.ztomo[p2]:
+ gamma_interp_z3 = 1
+ else:
+ gamma_interp_z3 = (z-self.ztomo[p2])/(self.ztomo[p6]-self.ztomo[p2])
+ if(gamma_interp_z3 > 1.): gamma_interp_z3 = 1.
+ if(gamma_interp_z3 < 0.): gamma_interp_z3 = 0.
+ #
+ if self.ztomo[p7]==self.ztomo[p3]:
+ gamma_interp_z4 = 1
+ else:
+ gamma_interp_z4 = (z-self.ztomo[p3])/(self.ztomo[p7]-self.ztomo[p3])
+ if(gamma_interp_z4 > 1.): gamma_interp_z4 = 1.
+ if(gamma_interp_z4 < 0.): gamma_interp_z4 = 0.
+ #
+ gamma_interp_z5 = 1. - gamma_interp_z1
+ gamma_interp_z6 = 1. - gamma_interp_z2
+ gamma_interp_z7 = 1. - gamma_interp_z3
+ gamma_interp_z8 = 1. - gamma_interp_z4
+ #
+ vp1 = self.vp[p0]
+ vp2 = self.vp[p1]
+ vp3 = self.vp[p2]
+ vp4 = self.vp[p3]
+ vp5 = self.vp[p4]
+ vp6 = self.vp[p5]
+ vp7 = self.vp[p6]
+ vp8 = self.vp[p7]
+ # [ ]
+ vs1 = self.vs[p0]
+ vs2 = self.vs[p1]
+ vs3 = self.vs[p2]
+ vs4 = self.vs[p3]
+ vs5 = self.vs[p4]
+ vs6 = self.vs[p5]
+ vs7 = self.vs[p6]
+ vs8 = self.vs[p7]
+ #
+ vp_final = vp1*(1.-gamma_interp_x)*(1.-gamma_interp_y)*(1.-gamma_interp_z1) + \
+ vp2*gamma_interp_x*(1.-gamma_interp_y)*(1.-gamma_interp_z2) + \
+ vp3*gamma_interp_x*gamma_interp_y*(1.-gamma_interp_z3) + \
+ vp4*(1.-gamma_interp_x)*gamma_interp_y*(1.-gamma_interp_z4) + \
+ vp5*(1.-gamma_interp_x)*(1.-gamma_interp_y)*gamma_interp_z1 + \
+ vp6*gamma_interp_x*(1.-gamma_interp_y)*gamma_interp_z2 + \
+ vp7*gamma_interp_x*gamma_interp_y*gamma_interp_z3 + \
+ vp8*(1.-gamma_interp_x)*gamma_interp_y*gamma_interp_z4
+ #
+ vs_final = vs1*(1.-gamma_interp_x)*(1.-gamma_interp_y)*(1.-gamma_interp_z1) + \
+ vs2*gamma_interp_x*(1.-gamma_interp_y)*(1.-gamma_interp_z2) + \
+ vs3*gamma_interp_x*gamma_interp_y*(1.-gamma_interp_z3) + \
+ vs4*(1.-gamma_interp_x)*gamma_interp_y*(1.-gamma_interp_z4) + \
+ vs5*(1.-gamma_interp_x)*(1.-gamma_interp_y)*gamma_interp_z1 + \
+ vs6*gamma_interp_x*(1.-gamma_interp_y)*gamma_interp_z2 + \
+ vs7*gamma_interp_x*gamma_interp_y*gamma_interp_z3 + \
+ vs8*(1.-gamma_interp_x)*gamma_interp_y*gamma_interp_z4
+ #
+ if(vp_final < self.vp_min): vp_final = self.vp_min
+ if(vs_final < self.vs_min): vs_final = self.vs_min
+ if(vp_final > self.vp_max): vp_final = self.vp_max
+ if(vs_final > self.vs_max): vs_final = self.vs_max
+ return vp_final,vs_final
+ else:
+ return vp_static,vs_static
+ #
+ def hex_simulation_parameter(self,h,vp_static=None,vs_static=None):
+ nodes=cubit.get_connectivity('Hex',h)
+ id_nodes=[0,1,2,3,4,5,6,7]
+ id_faces=[[1,3,4],[0,2,5],[1,3,6],[0,2,7],[5,7],[4,6],[5,7],[4,6]]
+ dt=[]
+ pmax=[]
+ vmax=[]
+ vmin=[]
+ for i in id_nodes[:-1]:
+ for j in id_nodes[i+1:]:
+ x1,y1,z1=cubit.get_nodal_coordinates(nodes[i])
+ x2,y2,z2=cubit.get_nodal_coordinates(nodes[j])
+ nvp1,nvs1=self.tomo(x1,y1,z1,vp_static,vs_static)
+ nvp2,nvs2=self.tomo(x2,y2,z2,vp_static,vs_static)
+ d=math.sqrt((x2-x1)**2+(y2-y1)**2+(z2-z1)**2)
+ #pmax_tmp=d*self.gllcoeff/min(nvp1,nvp2,nvs1,nvs2)*self.Ngll_per_wavelength
+ pmax_tmp=d*(.5-self.gllcoeff)/min(nvp1,nvp2,nvs1,nvs2)*self.Ngll_per_wavelength #more conservative.....
+ dt_tmp=self.Cmax*d*self.gllcoeff/max(nvp1,nvp2,nvs1,nvs2)
+ dt.append(dt_tmp)
+ pmax.append(pmax_tmp)
+ vmax.append(max(nvp1,nvp2,nvs1,nvs2))
+ vmin.append(min(nvp1,nvp2,nvs1,nvs2))
+ return min(dt),max(pmax),min(vmin),max(vmax)
+ #
+ def group_period(self):
+ tot=0.
+ if self.period_hyst is not None:
+ pmin=min(self.period_hyst.values())
+ period=self.hyst(pmin,self.period,self.period_hyst)
+ factor=(self.period-pmin)/self.nbin
+ for i in range(0,self.nbin+1):
+ if period.has_key(i):
+ txt='group "period_%.1e_%.1e" add hex ' %(pmin+factor*i,pmin+factor*(i+1))
+ txt=txt+' '.join(str(hh) for hh in period[i])
+ cubit.cmd(txt)
+ def group_timestep(self):
+ tot=0.
+ if self.dt_hyst is not None:
+ dtmax=max(self.dt_hyst.values())
+ dt=self.hyst(self.dt,dtmax,self.dt_hyst)
+ factor=(dtmax-self.dt)/self.nbin
+ for i in range(0,self.nbin+1):
+ if dt.has_key(i):
+ txt='group "timestep_%.1e_%.1e" add hex ' %(self.dt+factor*i,self.dt+factor*(i+1))
+ txt=txt+' '.join(str(hh) for hh in dt[i])
+ cubit.cmd(txt)
+
+
+#cubit.cmd('brick x 10000')
+#cubit.cmd('mesh vol 1')
+#cubit.cmd('refine hex in node in surf 1')
+#cubit.cmd('refine hex in node in surf 3')
+#vp=1000
+#vs=600
+#mesh=SEM_stability_3D()
+#mesh.check_simulation_parameter(vp_static=vp,vs_static=vs)
+#mesh.group_timestep()
+#mesh.group_period()
\ No newline at end of file
Added: seismo/3D/SPECFEM3D/trunk/CUBIT/local_volume.py
===================================================================
--- seismo/3D/SPECFEM3D/trunk/CUBIT/local_volume.py (rev 0)
+++ seismo/3D/SPECFEM3D/trunk/CUBIT/local_volume.py 2013-09-17 16:14:58 UTC (rev 22796)
@@ -0,0 +1,242 @@
+#############################################################################
+# local_volume.py
+# this file is part of GEOCUBIT #
+# #
+# Created by Emanuele Casarotti #
+# Copyright (c) 2008 Istituto Nazionale di Geofisica e Vulcanologia #
+# #
+#############################################################################
+# #
+# GEOCUBIT is free software: you can redistribute it and/or modify #
+# it under the terms of the GNU General Public License as published by #
+# the Free Software Foundation, either version 3 of the License, or #
+# (at your option) any later version. #
+# #
+# GEOCUBIT is distributed in the hope that it will be useful, #
+# but WITHOUT ANY WARRANTY; without even the implied warranty of #
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the #
+# GNU General Public License for more details. #
+# #
+# You should have received a copy of the GNU General Public License #
+# along with GEOCUBIT. If not, see <http://www.gnu.org/licenses/>. #
+# #
+#############################################################################
+try:
+ import start as start
+ cubit = start.start_cubit()
+except:
+ try:
+ import cubit
+ except:
+ print 'error importing cubit, check if cubit is installed'
+ pass
+
+def read_grid(filename=None):
+ import sys
+ import start as start
+ mpiflag,iproc,numproc,mpi = start.start_mpi()
+ #
+ numpy = start.start_numpy()
+ cfg = start.start_cfg(filename=filename)
+ from utilities import geo2utm
+
+ #
+ if cfg.nx and cfg.ny:
+ nx=cfg.nx
+ ny=cfg.ny
+ if cfg.nstep:
+ nx=min(cfg.nx,int(cfg.nx/cfg.nstep)+1)
+ ny=min(cfg.ny,int(cfg.ny/cfg.nstep)+1)
+ nstep=cfg.nstep
+ else:
+ nstep=1
+ else:
+ try:
+ xstep=cfg.step
+ ystep=cfg.step
+ except:
+ xstep=cfg.xstep
+ ystep=cfg.ystep
+ nx= int((cfg.longitude_max-cfg.longitude_min)/xstep)+1
+ ny= int((cfg.latitude_max-cfg.latitude_min)/ystep)+1
+ nstep=1
+ #
+ elev=numpy.zeros([nx,ny,cfg.nz],float)
+ coordx=numpy.zeros([nx,ny],float)
+ coordy=numpy.zeros([nx,ny],float)
+ #
+ if cfg.bottomflat:
+ elev[:,:,0] = cfg.depth_bottom
+ bottomsurface=1
+ else:
+ bottomsurface=0
+ #
+ for inz in range(bottomsurface,cfg.nz):
+ try:
+ grdfile = open(cfg.filename[inz-bottomsurface], 'r')
+ print 'reading ',cfg.filename[inz-bottomsurface]
+ except:
+ txt='error reading: '+ str( cfg.filename[inz-bottomsurface] )
+ raise NameError, txt
+ #
+ icoord=0
+ for iy in range(0,ny):
+ for ix in range(0,nx):
+ txt=grdfile.readline()
+ try:
+ if len(txt) != 0:
+ x,y,z=map(float,txt.split())
+ if iy%nstep == 0 and ix%nstep == 0:
+ icoord=icoord+1
+ x_current,y_current=geo2utm(x,y,cfg.unit)
+ jx=min(nx-1,ix/nstep)
+ jy=min(ny-1,iy/nstep)
+ coordx[jx,jy]=x_current
+ coordy[jx,jy]=y_current
+ elev[jx,jy,inz]=z
+ except:
+ print 'error reading point ',iy*cfg.nx+ix,txt, cfg.filename[inz-bottomsurface], ' proc ',iproc
+ raise NameError, 'error reading point'
+ #
+ if (nx)*(ny) != icoord:
+ if iproc == 0: print 'error in the surface file '+cfg.filename[inz-bottomsurface]
+ if iproc == 0: print 'x points ' +str(nx)+ ' y points ' +str(ny)+ ' tot points '+str((nx)*(ny))
+ if iproc == 0: print 'points read in '+cfg.filename[inz-bottomsurface]+': '+str(icoord)
+ raise NameError
+
+ #if iproc == 0: print 'end of reading grd ascii file '+cfg.filename[inz-bottomsurface]+' '+str(icoord)+ ' points'
+ grdfile.close()
+
+
+ return coordx,coordy,elev,nx,ny
+
+
+def extract_volume(xmin,ymin,xmax,ymax,coordx,coordy,elev,nx,ny,filename=None):
+ import sys
+ import start as start
+ #
+ mpiflag,iproc,numproc,mpi = start.start_mpi()
+ #
+ numpy = start.start_numpy()
+ cfg = start.start_cfg(filename=filename)
+
+ from utilities import geo2utm
+ #
+ rxstep=coordx[1,0]-coordx[0,0]
+ rystep=coordy[0,1]-coordy[0,0]
+
+ nxmin_cpu=min(0,int((x0-cfg.xmin)/rxstep)+1-10)
+ nymin_cpu=min(0,int((y0-cfg.ymin)/rxstep)+1-10)
+ nxmax_cpu=min(nx,int((x0-cfg.xmin)/rystep)+1+10)
+ nymax_cpu=min(ny,int((y0-cfg.ymin)/rystep)+1+10)
+ #
+ #
+ icurve=0
+ isurf=0
+ ivertex=0
+ #
+ #create vertex
+ last_surface=cubit.get_last_id('surface')
+ for inz in range(0,cfg.nz):
+ if cfg.bottomflat and inz == 0: #bottom layer
+
+ x_current,y_current=geo2utm(coordx[nxmin_cpu,nymin_cpu],coordy[nxmin_cpu,nymin_cpu],cfg.unit)
+ cubitcommand= 'create vertex '+ str( x_current )+ ' ' + str( y_current) +' '+ str( cfg.depth_bottom )
+ cubit.cmd(cubitcommand)
+ #
+ x_current,y_current=geo2utm(coordx[nxmin_cpu,nymax_cpu],coordy[nxmin_cpu,nymax_cpu],cfg.unit)
+ 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],cfg.unit)
+ 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],cfg.unit)
+ cubitcommand= 'create vertex '+ str( x_current )+ ' ' + str( y_current) +' '+ str( cfg.depth_bottom )
+ cubit.cmd(cubitcommand)
+ #
+ cubitcommand= 'create surface vertex 1 2 3 4'
+ cubit.cmd(cubitcommand)
+ #
+ isurf = isurf + 1
+
+ else:
+ vertex=[]
+
+ for iy in range(nymin_cpu,nymax_cpu+1):
+ ivx=0
+ for ix in range(nxmin_cpu,nxmax_cpu+1):
+ zvertex=elev[ix,iy,inz]
+ x_current,y_current=geo2utm(coordx[ix,iy],coordy[ix,iy],cfg.unit)
+ #
+ vertex.append(' Position '+ str( x_current ) +' '+ str( y_current )+' '+ str( zvertex ) )
+ #
+ print iproc, 'vertex created....'
+ n=max(nx,ny)
+ uline=[]
+ vline=[]
+ iv=0
+
+ cubit.cmd("set info off")
+ cubit.cmd("set echo off")
+ cubit.cmd("set journal off")
+
+ for iy in range(0,nymax_cpu-nymin_cpu+1):
+ positionx=''
+ for ix in range(0,nxmax_cpu-nxmin_cpu+1):
+ positionx=positionx+vertex[iv]
+ iv=iv+1
+ command='create curve spline '+positionx
+ cubit.cmd(command)
+ uline.append( cubit.get_last_id("curve") )
+ for ix in range(0,nxmax_cpu-nxmin_cpu+1):
+ positiony=''
+ for iy in range(0,nymax_cpu-nymin_cpu+1):
+ positiony=positiony+vertex[ix+iy*(nxmax_cpu-nxmin_cpu+1)]
+ command='create curve spline '+positiony
+ cubit.cmd(command)
+ vline.append( cubit.get_last_id("curve") )
+ #
+ cubit.cmd("set info "+cfg.cubit_info)
+ cubit.cmd("set echo "+cfg.echo_info)
+ cubit.cmd("set journal "+cfg.jou_info)
+ #
+ #
+ print iproc,'line created....'
+ umax=max(uline)
+ umin=min(uline)
+ vmax=max(vline)
+ vmin=min(vline)
+ cubitcommand= 'create surface net u curve '+ str( umin )+' to '+str( umax )+ ' v curve '+ str( vmin )+ ' to '+str( vmax )+' heal'
+ cubit.cmd(cubitcommand)
+ command = "del curve all"
+ cubit.cmd(command)
+ isurf=isurf+1
+ #
+ #
+ cubitcommand= 'del vertex all'
+ cubit.cmd(cubitcommand)
+ #cubit_error_stop(iproc,cubitcommand,ner)
+ cubitcommand= 'del curve all'
+ cubit.cmd(cubitcommand)
+ #
+ last_surface_2=cubit.get_last_id('surface')
+ #
+ for inz in range(1,cfg.nz):
+ #!cubit cmd
+ cubitcommand= 'create volume loft surface '+ str( inz+1 )+' '+str( inz )
+ cubit.cmd(cubitcommand)
+ #cubit_error_stop(iproc,cubitcommand,ner)
+ isurf=isurf+6
+ cubitcommand= 'del surface '+str(last_surface+1)+' to '+ str( last_surface_2 )
+ cubit.cmd(cubitcommand)
+ #cubit_error_stop(iproc,cubitcommand,ner)
+ #
+ #
+ cubit.cmd("set info "+cfg.cubit_info)
+ cubit.cmd("set echo "+cfg.echo_info)
+ cubit.cmd("set journal "+cfg.jou_info)
+ command = "compress all"
+ cubit.cmd(command)
+
\ No newline at end of file
Added: seismo/3D/SPECFEM3D/trunk/CUBIT/menu.py
===================================================================
--- seismo/3D/SPECFEM3D/trunk/CUBIT/menu.py (rev 0)
+++ seismo/3D/SPECFEM3D/trunk/CUBIT/menu.py 2013-09-17 16:14:58 UTC (rev 22796)
@@ -0,0 +1,315 @@
+#############################################################################
+# menu.py
+# this file is part of GEOCUBIT #
+# #
+# Created by Emanuele Casarotti #
+# Copyright (c) 2008 Istituto Nazionale di Geofisica e Vulcanologia #
+# #
+#############################################################################
+# #
+# GEOCUBIT is free software: you can redistribute it and/or modify #
+# it under the terms of the GNU General Public License as published by #
+# the Free Software Foundation, either version 3 of the License, or #
+# (at your option) any later version. #
+# #
+# GEOCUBIT is distributed in the hope that it will be useful, #
+# but WITHOUT ANY WARRANTY; without even the implied warranty of #
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the #
+# GNU General Public License for more details. #
+# #
+# You should have received a copy of the GNU General Public License #
+# along with GEOCUBIT. If not, see <http://www.gnu.org/licenses/>. #
+# #
+#############################################################################
+
+import getopt, sys
+
+def usage():
+ print """
+ GEOCUBIT HELP...
+
+ 1) UTILITIES
+
+ check the configuration of the libraries and dependencies:
+ GEOCUBIT.py --chklib
+
+ check the parameter file:
+ GEOCUBIT.py --chkcfg --cfg=[filename]
+
+ 2) CREATE GEOMETRY
+
+ create a surface from regular ascii grid or ascii lines defining a skin:
+ GEOCUBIT.py --surface=[surface file] (--regulargrid=[options]) (--skin=[options])
+
+ create a plane surface
+ GEOCUBIT.py --plane --x1=[x,y,z] --x2=[x,y,z] --x3=[x,y,z] --x4=[x,y,z] --unit=[utm/geo]
+
+ create acis surfaces using a parameter file:
+ GEOCUBIT.py --build_surface --cfg=[filename]
+
+ **SERIAL**: create cubit volumes using a parameter file:
+ GEOCUBIT.py --build_volume --cfg=[filename] (--id_proc=[num_processor, default=0])
+
+ **PARALLEL**: create a volume from a parameter file
+ mpirun -n [numproc] pyMPI GEOCUBIT.py --build_volume --cfg=[filename
+
+ 3) MESHING
+
+ **SERIAL**: meshing a volumes
+ GEOCUBIT.py --mesh --cfg=[filename] (--id_proc=[num_processor, default=0]) - note: without the --build_volume flag the script recall an old 'geometry_vol_[id_proc].cub' file
+
+ build a volume and mesh it....
+ GEOCUBIT.py --build_volume --mesh --cfg=[filename] (--id_proc=[num_processor, default=0])
+
+ **PARALLEL**: meshing a volume from a parameter file (it is possible to skip the build_volume and use some geometry files already created, but only with a parallel file system and it is risky)
+ mpirun -n [numproc] pyMPI GEOCUBIT.py --build_volume --mesh --cfg=[filename]
+
+ 4) FINALIZING AND EXPORTING
+
+ collect some cubit files and merge in a single free mesh cubitfile
+ GEOCUBIT.py --collect --merge --meshfiles=[list of files] --cpux=N --cpuy=N (--rangecpux=[cpuxmin,cpuxmax], --rangecpuy=[cpuymin,cpuymax])
+
+ collect a single free mesh cubitfile and refine the hex inside some curve (ex. basin)
+ GEOCUBIT.py --collect --meshfiles=[list of files] --curverefining=[list of SAT files]
+
+ export a cubit mesh file (with blocks defined following the note) in a SPECFEM3D_SESAME mesh
+ GEOCUBIT.py --export2SPECFEM3D --meshfiles=[filename] (--listblock=block1,block2,..,blockN --listflag=[list of specfem flag, i.e. --listflag=1,2,3,-1])
+
+ """
+try:
+ opts, args = getopt.getopt(sys.argv[1:], "sjmohbp1", ["SEMoutput=","qlog","mfast","curverefining=","output=","rangecpux=","rangecpuy=","equivalence","listflag=","listblock=","cpux=","cpuy=","exofiles=","partitioner","plane","x1=","x2=","x3=","x4=","unit=","chkcfg","mat=","merge_tolerance=","export2SPECFEM3D=","mesh","chklib","cfg=","job=","basin","help", "id_proc=", "surface=","script","jou","strat","MPI","regulargrid=",'skin=',"build_surface","build_volume","merge1","merge2","merge","collect","meshfiles="])
+ print opts, args
+except getopt.GetoptError,errmsg:
+ if str(errmsg) == 'option --export2SPECFEM3D requires argument':
+ for i,xop in enumerate(sys.argv[1:]):
+ if 'export2SPECFEM3D' in xop:
+ sys.argv[i+1]=sys.argv[i+1]+'=.'
+ try:
+ opts, args = getopt.getopt(sys.argv[1:], "sjmohbp1", ["SEMoutput=","qlog","mfast","curverefining=","output=","rangecpux=","rangecpuy=","equivalence","listflag=","listblock=","cpux=","cpuy=","exofiles=","partitioner","plane","x1=","x2=","x3=","x4=","unit=","chkcfg","mat=","merge_tolerance=","export2SPECFEM3D=","mesh","chklib","cfg=","job=","basin","help", "id_proc=", "surface=","script","jou","strat","MPI","regulargrid=",'skin=',"build_surface","build_volume","merge1","merge2","merge","collect","meshfiles="])
+ print opts, args
+ except getopt.GetoptError,errmsg:
+ print str(errmsg)
+ usage()
+ sys.exit(2)
+ else:
+ print str(errmsg)
+ usage()
+ sys.exit()
+except AttributeError:
+ opts=None
+ args=None
+ print opts, args
+
+output='totalmesh_merged'
+SPECFEM3D_output_dir='.'
+verbose = False
+surface = False
+script = False
+mesh = False
+basin = False
+configuration = None
+id_proc=0
+single=False
+nomesh=False
+jobid=0
+build_surface=False
+build_volume=False
+meshing=False
+ckbound_method1=False
+ckbound_method2=False
+collect=False
+export2SPECFEM3D=False
+merge_tolerance=0
+material_file='material_archive.dat'
+material_assignement=[]
+chkcfg=False
+create_plane=False
+create_partitioner=False
+cubfiles=None
+exofiles=None
+listflag=None
+listblock=None
+cpux=0
+cpuy=0
+cpuxmin=0
+cpuymin=0
+cpuxmax=None
+cpuymax=None
+curverefining=False
+
+
+
+
+
+qlog=False
+
+
+if opts:
+ for o, value in opts:
+ #print o,value
+ if o in ('--partitioner'):
+ create_partitioner=True
+ if o == ('--surface'):
+ surface=True
+ surface_name=value
+ if o == ('--build_surface'):
+ build_surface=True
+ if o == ('--build_volume'):
+ build_volume=True
+ if surface and o == ('--regular_grid'):
+ surface_type='regular_grid'
+ tmp=value.split('/')
+ if len(tmp)==4:
+ num_x,num_y,unit,delimiter=value.split('/')
+ elif len(tmp)==3:
+ num_x,num_y,unit=value.split('/')
+ num_x=int(num_x)
+ num_y=int(num_y)
+ delimiter=' '
+ if surface and o == ('--skin'):
+ surface_type='skin'
+ tmp=value.split('/')
+ if len(tmp)==4:
+ directionx,directiony,unit,delimiter=value.split('/')
+ elif len(tmp)==3:
+ directionx,directiony,unit=value.split('/')
+ directiony=int(directiony)
+ directionx=int(directionx)
+ if o in ('--plane'):
+ create_plane=True
+ cfg_name=False
+ if o in ('--x1'):
+ x1=value
+ if o in ('--x2'):
+ x2=value
+ if o in ('--x3'):
+ x3=value
+ if o in ('--x4'):
+ x4=value
+ if o in ('--unit'):
+ unit=value
+ #
+ if o in ('--build_volume'):
+ build_volume=True
+ if o in ("-h", "--help"):
+ usage()
+ sys.exit(2)
+ #chec k the configuration
+ if o in ("--chklib"):
+ import start as start
+ mpiflag,iproc,numproc,mpi = start.start_mpi()
+ if mpiflag:
+ print '--------, MPI ON, parallel mesher ready'
+ else:
+ print '--------, MPI OFF, serial mesher ready'
+ numpy = start.start_numpy()
+ print '--------, Numpy ON'
+ cubit = start.start_cubit()
+ print '--------, CUBIT ON'
+ sys.exit()
+ if o in ("--cfg"):
+ cfg_name=value
+ if o == ('--surface'):
+ surface=True
+ surface_name=value
+ if o == ("--mesh"):
+ meshing=True
+ if o in ("-1"):
+ single=True
+ if o in ("--collect"):
+ collect=True
+ if o in ("--merge2") and o != '--merge':
+ ckbound_method2=True
+ if o in ("--equivalence","--merge1","--merge"):
+ ckbound_method1=True
+ if o in ("--mfast"):
+ ckbound_method1=True
+ ckbound_method2=True
+ if o in ("--meshfiles"):
+ cubfiles=value
+ if o in ("--exofiles"):
+ exofiles=value
+ if o in ("--export2SPECFEM3D"):
+ export2SPECFEM3D=True
+ SPECFEM3D_output_dir=value
+ import os
+ try:
+ os.makedirs(SPECFEM3D_output_dir)
+ except OSError:
+ pass
+ if o in ("--merge_tolerance") and o != '--merge' and o != '--merge2' and o != '--merge1':
+ merge_tolerance=map(float,value.split(','))
+ if o in ("--mat"):
+ material_assignement.append([value.split(',')[0],value.split(',')[1]])
+ if o in ("--listblock"):
+ listblock=map(int,value.split(','))
+ if o in ("--listflag"):
+ listflag=map(int,value.split(','))
+ if o in ("--chkcfg"):
+ chkcfg=True
+ if o in ('--id_proc'):
+ id_proc=int(value)
+ if o in ('--rangecpux'):
+ cpuxmin=int(value.split(',')[0])
+ cpuxmax=int(value.split(',')[1])+1
+ if o in ('--rangecpuy'):
+ cpuymin=int(value.split(',')[0])
+ cpuymax=int(value.split(',')[1])+1
+ if o in ('--cpux'):
+ cpux=int(value)
+ if o in ('--cpuy'):
+ cpuy=int(value)
+ if o in ("--qlog"):
+ qlog=True
+ if o in ("--output","-o"):
+ output=value
+ if o in ("--curverefining"):
+ curverefining=value.split(',')
+ if o in ("SEMoutput"):
+ SPECFEM3D_output_dir=value
+ print cpuxmax,cpuymax
+ if cpuymax:
+ pass
+ elif cpuy > 1:
+ cpuymax=cpuy
+ else:
+ cpuymax=1
+ if cpuxmax:
+ pass
+ elif cpux > 1:
+ cpuxmax=cpux
+ else:
+ cpuxmax=1
+ print cpuxmax,cpuymax
+ if chkcfg==True:
+ import start as start
+ cfg=start.start_cfg()
+ d=cfg.__dict__
+ ks=d.keys()
+ ks.sort()
+ for k in ks:
+ if '__' not in k and '<' not in str(d[k]) and d[k] is not None:
+ txt=str(k)+' -----> '+str(d[k]) #
+ txt=txt.replace("'","").replace('"','') #
+ print txt #
+ else:
+ try:
+ import start as start
+ cfg=start.start_cfg()
+ f=open('cfg.log','w')
+ print>>f, 'CFG FILE: ',cfg_name
+ d=cfg.__dict__
+ ks=d.keys()
+ ks.sort()
+ for k in ks:
+ if '__' not in k and '<' not in str(d[k]) and d[k] is not None:
+ txt=str(k)+' -----> '+str(d[k]) #
+ txt=txt.replace("'","").replace('"','') #
+ print>>f, txt
+ f.close() #
+ except:
+ pass
+elif opts == []:
+ print __name__
+ usage()
+ raise AttributeError('no options')
Added: seismo/3D/SPECFEM3D/trunk/CUBIT/mesh_volume.py
===================================================================
--- seismo/3D/SPECFEM3D/trunk/CUBIT/mesh_volume.py (rev 0)
+++ seismo/3D/SPECFEM3D/trunk/CUBIT/mesh_volume.py 2013-09-17 16:14:58 UTC (rev 22796)
@@ -0,0 +1,443 @@
+#############################################################################
+# mesh_volume.py
+# this file is part of GEOCUBIT #
+# #
+# Created by Emanuele Casarotti #
+# Copyright (c) 2008 Istituto Nazionale di Geofisica e Vulcanologia #
+# #
+#############################################################################
+# #
+# GEOCUBIT is free software: you can redistribute it and/or modify #
+# it under the terms of the GNU General Public License as published by #
+# the Free Software Foundation, either version 3 of the License, or #
+# (at your option) any later version. #
+# #
+# GEOCUBIT is distributed in the hope that it will be useful, #
+# but WITHOUT ANY WARRANTY; without even the implied warranty of #
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the #
+# GNU General Public License for more details. #
+# #
+# You should have received a copy of the GNU General Public License #
+# along with GEOCUBIT. If not, see <http://www.gnu.org/licenses/>. #
+# #
+#############################################################################
+
+try:
+ import start as start
+ cubit = start.start_cubit()
+except:
+ try:
+ import cubit
+ except:
+ print 'error importing cubit, check if cubit is installed'
+ pass
+
+
+def mesh(filename=None):
+ """create the mesh"""
+ import start as start
+ cfg = start.start_cfg(filename=filename)
+ mpiflag,iproc,numproc,mpi = start.start_mpi()
+ #
+ if cfg.map_meshing_type == 'regularmap':
+ mesh_layercake_regularmap(filename=filename)
+ #elif cfg.map_meshing_type == 'partitioner':
+ # mesh_partitioner()
+ else:
+ print 'error: map_meshing_type ', cfg.map_meshing_type,' not implemented'
+
+
+def mesh_layercake_regularmap(filename=None):
+ import sys,os
+ import start as start
+ #
+ class cubitvolume:
+ def __init__(self,ID,intervalv,centerpoint,dimension):
+ self.ID=ID
+ self.intervalv=intervalv
+ self.centerpoint=centerpoint
+ self.dim=dimension
+
+ def __repr__(self):
+ msg="(vol:%3i, vertical interval: %4i, centerpoint: %8.2f)" % (self.ID, self.intervalv,self.centerpoint)
+ return msg
+ #
+ def by_z(x,y):
+ return cmp(x.centerpoint,y.centerpoint)
+ #
+ #
+ mpiflag,iproc,numproc,mpi = start.start_mpi()
+ from utilities import importgeometry,savemesh,get_v_h_list
+ #
+ numpy = start.start_numpy()
+ cfg = start.start_cfg(filename=filename)
+ from math import sqrt
+ from sets import Set
+ #
+ list_vol=cubit.parse_cubit_list("volume","all")
+ if len(list_vol) != 0:
+ pass
+ else:
+ geometryfile='geometry_vol_'+str(iproc)+'.cub'
+ importgeometry(geometryfile,iproc=iproc)
+ #
+ command = 'composite create curve all'
+ cubit.cmd(command)
+ print 'NO CRITICAL ERROR: "No valid composites can be created from the specified curves." is not critical. \n It means that your model is clean and you don"t need a virtual geometry'
+ #
+ command = "compress all"
+ cubit.cmd(command)
+ list_vol=cubit.parse_cubit_list("volume","all")
+ nvol=len(list_vol)
+ vol=[]
+ for id_vol in list_vol:
+ p=cubit.get_center_point("volume",id_vol)
+ vol.append(cubitvolume(id_vol,1,p[2],0))
+ vol.sort(by_z)
+ #
+ for id_vol in range(0,nvol):
+ vol[id_vol].intervalv=cfg.iv_interval[id_vol]
+ #
+ #
+ surf_vertical=[]
+ surf_or=[]
+ top_surface=0
+ top_surface_add=''
+ bottom_surface=0
+ #
+ zmin_box=cubit.get_total_bounding_box("volume",list_vol)[6]
+ xmin_box=cubit.get_total_bounding_box("volume",list_vol)[0]
+ xmax_box=cubit.get_total_bounding_box("volume",list_vol)[1]
+ ymin_box=cubit.get_total_bounding_box("volume",list_vol)[3]
+ ymax_box=cubit.get_total_bounding_box("volume",list_vol)[4]
+ #
+ #
+ #interval assignement
+ surf_or,surf_vertical,list_curve_or,list_curve_vertical,bottom,top = get_v_h_list(list_vol)
+
+ for k in surf_vertical:
+ command = "surface "+str(k)+" scheme submap"
+ cubit.cmd(command)
+ for k in surf_or:
+ command = "surface "+str(k)+" scheme "+cfg.or_mesh_scheme
+ cubit.cmd(command)
+ for k in list_curve_or:
+ length=cubit.get_curve_length(k)
+ interval=int(2*round(.5*length/cfg.size,0))
+ command = "curve "+str(k)+" interval "+str(interval)
+ cubit.cmd(command)
+ #cubit_error_stop(iproc,command,ner)
+ command = "curve "+str(k)+" scheme equal"
+ cubit.cmd(command)
+ #cubit_error_stop(iproc,command,ner)
+ #
+ for s in surf_vertical:
+ lcurve=cubit.get_relatives("surface",s,"curve")
+ interval_store=[]
+ for k in lcurve:
+ interval_curve=cubit.get_mesh_intervals('curve',k)
+ if k in list_curve_vertical:
+ volume_id = cubit.get_owning_volume("curve", k)
+ for idv in range(0,nvol):
+ if vol[idv].ID == volume_id:
+ int_v=vol[idv].intervalv
+ command = "curve "+str(k)+" interval "+str(int_v)
+ cubit.cmd(command)
+ #cubit_error_stop(iproc,command,ner)
+ command = "curve "+str(k)+" scheme equal"
+ cubit.cmd(command)
+ #cubit_error_stop(iproc,command,ner)
+ else:
+ interval_store.append((k,interval_curve))
+ if len(interval_store) != 0:
+ interval_min=min([iv[1] for iv in interval_store])
+ command = "curve "+' '.join(str(iv[0]) for iv in interval_store)+" interval "+str(interval_min)
+ cubit.cmd(command)
+ #cubit_error_stop(iproc,command,ner)
+ command = "curve "+' '.join(str(iv[0]) for iv in interval_store)+" scheme equal"
+ cubit.cmd(command)
+ #cubit_error_stop(iproc,command,ner)
+ #cubit_error_stop(iproc,command,ner)
+ #
+ #meshing
+ if cfg.or_mesh_scheme == 'pave':
+ command='mesh surf '+' '.join(str(t) for t in top)
+ cubit.cmd(command)
+ elif cfg.or_mesh_scheme == 'map':
+ command='mesh surf '+' '.join(str(t) for t in bottom)
+ cubit.cmd(command)
+ for id_volume in range(nvol-1,-1,-1):
+ command = "mesh vol "+str(vol[id_volume].ID)
+ cubit.cmd(command)
+
+ #
+ #smoothing
+ 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)
+ cubitcommand= 'surf all smooth scheme laplacian '
+ cubit.cmd(cubitcommand)
+ cubitcommand= 'smooth surf all'
+ cubit.cmd(cubitcommand)
+ #
+ cubitcommand= 'vol all smooth scheme laplacian '
+ cubit.cmd(cubitcommand)
+ cubitcommand= 'smooth vol all'
+ cubit.cmd(cubitcommand)
+ #
+ #
+ ##vertical refinement
+ ##for nvol = 3
+ ##
+ ##___________________________ interface 4
+ ##
+ ##vol 2
+ ##___________________________ interface 3
+ ##
+ ##vol 1
+ ##___________________________ interface 2
+ ##
+ ##vol 0
+ ##___________________________ interface 1
+ ##
+ refinement(nvol,vol,filename=filename)
+ #
+ #top layer vertical coarsening
+ print 'coarsening top layer... ',cfg.coarsening_top_layer
+ if cfg.coarsening_top_layer:
+ from sets import Set
+ cubitcommand= 'del mesh vol '+str(vol[-1].ID)+ ' propagate'
+ cubit.cmd(cubitcommand)
+ s1=Set(list_curve_vertical)
+ print s1
+ command = "group 'list_curve_tmp' add curve "+"in vol "+str(vol[-1].ID)
+ cubit.cmd(command)
+ group=cubit.get_id_from_name("list_curve_tmp")
+ list_curve_tmp=cubit.get_group_curves(group)
+ command = "delete group "+ str(group)
+ cubit.cmd(command)
+ s2=Set(list_curve_tmp)
+ print s2
+ lc=list(s1 & s2)
+ print lc
+ #
+ cubitcommand= 'curve '+' '.join(str(x) for x in lc)+' interval '+str(cfg.actual_vertical_interval_top_layer)
+ cubit.cmd(cubitcommand)
+ cubitcommand= 'mesh vol '+str(vol[-1].ID)
+ cubit.cmd(cubitcommand)
+ #
+ n=cubit.get_sideset_id_list()
+ if len(n) != 0:
+ command = "del sideset all"
+ cubit.cmd(command)
+ n=cubit.get_block_id_list()
+ if len(n) != 0:
+ command = "del block all"
+ cubit.cmd(command)
+ #
+ import boundary_definition
+ entities=['face']
+ print iproc, 'hex block definition...'
+ boundary_definition.define_bc(entities,parallel=True,cpux=cfg.cpux,cpuy=cfg.cpuy,cpuxmin=0,cpuymin=0)
+ #save mesh
+
+ 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...'
+ savemesh(mpiflag,iproc=iproc,filename=filename)
+ #
+
+def refinement(nvol,vol,filename=None):
+ import start as start
+ cfg = start.start_cfg(filename=filename)
+ #
+ #vertical refinement
+ #for nvol = 3
+ #
+ #___________________________ interface 4
+ #
+ #vol 2
+ #___________________________ interface 3
+ #
+ #vol 1
+ #___________________________ interface 2
+ #
+ #vol 0
+ #___________________________ interface 1
+ #
+ #
+ if cfg.ntripl != 0:
+ if len(cfg.refinement_depth) != 0:
+ #get the topo surface....
+ surf=cubit.get_relatives('volume',vol[nvol+1-2].ID,'surface')
+ zstore=[-1,-999999999]
+ for s in surf:
+ c=cubit.get_center_point('surface',s)
+ z=c[2]
+ if z > zstore[1]:
+ zstore=[s,z]
+ tsurf=zstore[0]
+ for idepth in cfg.refinement_depth:
+ cubitcommand= 'refine node in surf '+str(tsurf)+' numsplit 1 bias 1.0 depth '+str(idepth)
+ cubit.cmd(cubitcommand)
+ else:
+ for ir in cfg.tripl:
+ if ir == 1:
+ command = "comment '"+"interface = 1 means that the refinement interface is at the bottom of the volume"+"'"
+ cubit.cmd(command)
+ txt=' all '
+ idepth = 1
+ cubitcommand= 'refine hex in vol '+txt
+ elif ir != nvol+1:
+ txt=''
+ for id_vol_ref in range(ir-1,nvol):
+ txt=txt+str(vol[id_vol_ref].ID)+' '
+ #txt=txt+'except hex in vol '+str(vol[ir-2].ID)
+ #idepth = 1
+ #try:
+ # if cfg.refine_basin:
+ # idepth=2
+ #except:
+ # pass
+ cubitcommand= 'refine hex in vol '+txt
+ else:
+ #refinement on the top surface
+
+ surf=cubit.get_relatives('volume',vol[ir-2].ID,'surface')
+ zstore=[-1,-999999999]
+ for s in surf:
+ c=cubit.get_center_point('surface',s)
+ z=c[2]
+ if z > zstore[1]:
+ zstore=[s,z]
+ idepth=1
+ cubitcommand= 'refine node in surf '+str(zstore[0])+' numsplit 1 bias 1.0 depth '+str(idepth)
+ cubit.cmd(cubitcommand)
+
+
+
+
+def pinpoly(x,y,polyx,polyy):
+ """point in polygon using ray tracing"""
+ n = len(polyx)
+ inside = False
+ poly=zip(polyx,polyy)
+ px_0,py_0 = poly[0]
+ for i in range(n+1):
+ px_1,py_1 = poly[i % n]
+ if y > min(py_0,py_1):
+ if y <= max(py_1,py_0):
+ if x <= max(px_1,px_0):
+ if py_0 != py_1:
+ intersect = (y-py_0)*(px_1-px_0)/(py_1-py_0)+px_0
+ if px_1 == px_0 or x <= intersect:
+ inside = not inside
+ px_0,py_0 = px_1,py_1
+ return inside
+
+def curve2poly(line):
+ curve=int(line)
+ cubit.cmd('curve '+str(curve)+' size auto factor 1')
+ cubit.cmd('mesh curve '+str(curve))
+ n=cubit.get_curve_nodes(curve)
+ orientnode=[]
+ vertex_list = cubit.get_relatives("curve", curve, "vertex")
+ if len(vertex_list) != 0:
+ startnode=cubit.get_vertex_node(vertex_list[0])
+ cubit.cmd('del group pgon')
+ cubit.cmd("group 'pgon' add edge in node "+str(startnode))
+ group1 = cubit.get_id_from_name("pgon")
+ edges = list(cubit.get_group_edges(group1))
+ edgestart=edges[0]
+ else:
+ startnode=n[0]
+ cubit.cmd('del group pgon')
+ cubit.cmd("group 'pgon' add edge in node "+str(startnode))
+ group1 = cubit.get_id_from_name("pgon")
+ edges = list(cubit.get_group_edges(group1))
+ edgestart=edges[0]
+ begin=startnode
+ orientnode.append(begin)
+ node_id_list = list(cubit.get_connectivity("edge", edgestart))
+ node_id_list.remove(startnode)
+ startnode=node_id_list[0]
+ orientnode.append(startnode)
+ stopflag=False
+ while startnode != begin and not stopflag:
+ cubit.cmd('del group pgon')
+ cubit.cmd("group 'pgon' add edge in node "+str(startnode))
+ group1 = cubit.get_id_from_name("pgon")
+ edges = list(cubit.get_group_edges(group1))
+ if len(edges) != 1:
+ edges.remove(edgestart)
+ edgestart=edges[0]
+ node_id_list = list(cubit.get_connectivity("edge", edgestart))
+ node_id_list.remove(startnode)
+ orientnode.append(node_id_list[0])
+ startnode=node_id_list[0]
+ else:
+ stopflag=True
+ vx=[]
+ vy=[]
+ for n in orientnode:
+ v=cubit.get_nodal_coordinates(n)
+ vx.append(v[0])
+ vy.append(v[1])
+ return vx,vy,orientnode
+
+def get_nodes_inside_curve(nodes, curve):
+ vx,vy,n=curve2poly(curve)
+ nodes_inside=[]
+ for n in nodes:
+ vp=cubit.get_nodal_coordinates(n)
+ if pinpoly(vp[0],vp[1],vx,vy): nodes_inside.append(n)
+ return nodes_inside
+
+def refine_inside_curve(curves,ntimes=1,depth=1,block=1,surface=False):
+ if not isinstance(curves,list):
+ if isinstance(curves,str):
+ curves=map(int,curves.split())
+ else:
+ curves=[curves]
+ for curve in curves:
+ cubit.cmd('del group ntop')
+ if not surface:
+ cubit.cmd("group 'ntop' add node in face in block "+str(block))
+ else:
+ cubit.cmd("group 'ntop' add node in face in surface "+str(surface))
+ group1 = cubit.get_id_from_name("ntop")
+ nodes = list(cubit.get_group_nodes(group1))
+ ni=get_nodes_inside_curve(nodes, curve)
+ if ntimes > 1:
+ cmd='del group hex_refining'
+ cubit.cmd(cmd)
+ command = "group 'hex_refining' add hex propagate face in node "+' '.join(str(x) for x in ni)+"times "+str(ntimes)
+ cubit.cmd(command)
+ id_group=cubit.get_id_from_name('hex_refining')
+ command='refine hex in group "hex_refining" numsplit 1 bias 1.0 depth '+str(depth)+' smooth'
+ cubit.cmd(command)
+ else:
+ command='refine node '+" ".join(str(x) for x in ni)+' numsplit 1 bias 1.0 depth '+str(depth)+' smooth'
+ cubit.cmd(command)
+ #
+ #
+ cmd='group "negativejac" add quality hex all Jacobian high'
+ cubit.cmd(cmd)
+ 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'
+ import sys
+ #sys.exit()
+
Added: seismo/3D/SPECFEM3D/trunk/CUBIT/quality_log.py
===================================================================
--- seismo/3D/SPECFEM3D/trunk/CUBIT/quality_log.py (rev 0)
+++ seismo/3D/SPECFEM3D/trunk/CUBIT/quality_log.py 2013-09-17 16:14:58 UTC (rev 22796)
@@ -0,0 +1,282 @@
+#############################################################################
+# quality_log.py
+# this file is part of GEOCUBIT #
+# #
+# Created by Emanuele Casarotti #
+# Copyright (c) 2008 Istituto Nazionale di Geofisica e Vulcanologia #
+# #
+#############################################################################
+# #
+# GEOCUBIT is free software: you can redistribute it and/or modify #
+# it under the terms of the GNU General Public License as published by #
+# the Free Software Foundation, either version 3 of the License, or #
+# (at your option) any later version. #
+# #
+# GEOCUBIT is distributed in the hope that it will be useful, #
+# but WITHOUT ANY WARRANTY; without even the implied warranty of #
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the #
+# GNU General Public License for more details. #
+# #
+# You should have received a copy of the GNU General Public License #
+# along with GEOCUBIT. If not, see <http://www.gnu.org/licenses/>. #
+# #
+#############################################################################
+try:
+ import start as start
+ cubit = start.start_cubit()
+except:
+ try:
+ import cubit
+ except:
+ print 'error importing cubit, check if cubit is installed'
+ pass
+
+
+def quality_log(tqfile=None):
+ """
+ creation of the quality parameter file
+ """
+ import start as start
+ #
+ #
+ mpiflag,iproc,numproc,mpi = start.start_mpi()
+ #
+ #
+ from hex_metric import SEM_metric_3D
+ #
+ lvol=cubit.parse_cubit_list('volume','all')
+ if len(lvol)!=0:
+ cubit.cmd('quality vol all allmetric ')
+ else:
+ cubit.cmd('quality hex in block all allmetric ')
+ cubit.cmd('list model ')
+ #
+ toclose=True
+ if isinstance(tqfile,file):
+ totstat_file=tqfile
+ elif isinstance(tqfile,str):
+ totstat_file=open(tqfile+'_cubitquality_skewness_proc_'+str(iproc)+'.log','w')
+ else:
+ import sys
+ totstat_file=sys.stdout
+ toclose=False
+
+ mesh=SEM_metric_3D()
+ mesh.check_metric()
+
+ if mesh.max_skewness is not None:
+ mesh.skew_hystogram=mesh.hyst(0,mesh.max_skewness,mesh.skew_hyst)
+ totstat_file.write('-'*70+'\n')
+
+
+
+
+
+ totstat_file.write('='*70+'\n')
+ totstat_file.write('SKEWNESS'+'\n')
+ totstat_file.write('='*70+'\n')
+ if len(mesh.hex_max_skewness) <= 30:
+ totstat_file.write('max = '+str(mesh.max_skewness)+' in hexes '+str(mesh.hex_max_skewness)+'\n')
+ totstat_file.write('(angle -> minimun ='+str(mesh.min_angle)+ ' maximun ='+str(mesh.max_angle)+')'+'\n')
+ else:
+ totstat_file.write('max = '+str(mesh.max_skewness)+' in '+str(len(mesh.hex_max_skewness))+' hexes '+'\n')
+ totstat_file.write('(angle -> minimun ='+str(mesh.min_angle)+' maximun ='+str(mesh.max_angle)+')'+'\n')
+ totstat_file.write('-'*70+'\n')
+ totstat_file.write('skew hystogram')
+ totstat_file.write('-'*70+'\n')
+ tot=0
+ for i in mesh.skew_hystogram.values():
+ tot=tot+len(i)
+ #k=mesh.skew_hystogram.keys()
+ #k.sort()
+ factor=mesh.max_skewness/mesh.nbin
+ for i in range(0,mesh.nbin+1):
+ if mesh.skew_hystogram.has_key(i):
+ if (i+1)*factor <= 1:
+ totstat_file.write(str(i)+' ['+str(i*factor)+'->'+str((i+1)*factor)+'[ : '+str(len(mesh.skew_hystogram[i]))+'/'+str(tot)+' hexes ('+str(len(mesh.skew_hystogram[i])/float(tot)*100.)+'%)'+'\n')
+ else:
+ if (i+1)*factor <= 1:
+ totstat_file.write(str(i)+' ['+str(i*factor)+'->'+str((i+1)*factor)+'[ : 0/'+str(tot)+' hexes (0%)'+'\n')
+ totstat_file.write('-'*70+'\n')
+ ###############################################
+ if mesh.min_edge_length is not None:
+ mesh.edgemin_hystogram=mesh.hyst(mesh.min_edge_length,mesh.max_edge_length,mesh.edgemin_hyst)
+ mesh.edgemax_hystogram=mesh.hyst(mesh.min_edge_length,mesh.max_edge_length,mesh.edgemax_hyst)
+ totstat_file.write('='*70+'\n')
+ totstat_file.write('edge length')
+ totstat_file.write('='*70+'\n')
+ if len(mesh.hex_min_edge_length) <= 30:
+ totstat_file.write('minimum edge length: '+str(mesh.min_edge_length)+ ' in hexes '+str(mesh.hex_min_edge_length)+'\n')
+ else:
+ totstat_file.write('minimum edge length: '+str(mesh.min_edge_length)+ ' in '+str(len(mesh.hex_min_edge_length))+ ' hexes.'+'\n')
+ if len(mesh.hex_max_edge_length) <= 30:
+ totstat_file.write('maximum edge length: '+str(mesh.max_edge_length)+ ' in hexes '+str(mesh.hex_max_edge_length)+'\n')
+ else:
+ totstat_file.write('maximum edge length: '+str(mesh.max_edge_length)+' in '+str(len(mesh.hex_max_edge_length))+ ' hexes.'+'\n')
+ totstat_file.write('-'*70+'\n')
+ totstat_file.write('edge length hystogram')
+ totstat_file.write('-'*70+'\n')
+ factor=(mesh.max_edge_length-mesh.min_edge_length)/mesh.nbin
+ totstat_file.write('minimum edge length'+'\n')
+ tot=0
+ for i in mesh.edgemin_hystogram.values():
+ tot=tot+len(i)
+ #k=mesh.edgemin_hystogram.keys()
+ #k.sort()
+ for i in range(0,mesh.nbin+1):
+ if mesh.edgemin_hystogram.has_key(i):
+ totstat_file.write(str(i)+' ['+str(i*factor+mesh.min_edge_length)+'->'+str((i+1)*factor+mesh.min_edge_length)+'[ : '+str(len(mesh.edgemin_hystogram[i]))+'/'+str(tot)+' hexes ('+str(len(mesh.edgemin_hystogram[i])/float(tot)*100.)+'%)'+'\n')
+ else:
+ totstat_file.write(str(i)+' ['+str(i*factor+mesh.min_edge_length)+'->'+str((i+1)*factor+mesh.min_edge_length)+'[ : 0/'+str(tot)+' hexes (0%)'+'\n')
+ totstat_file.write('-'*70+'\n')
+ totstat_file.write('maximum edge length')
+ tot=0
+ for i in mesh.edgemax_hystogram.values():
+ tot=tot+len(i)
+ #k=mesh.edgemax_hystogram.keys()
+ #k.sort()
+ for i in range(0,mesh.nbin+1):
+ if mesh.edgemax_hystogram.has_key(i):
+ totstat_file.write(str(i)+' ['+str(i*factor+mesh.min_edge_length)+'->'+str((i+1)*factor+mesh.min_edge_length)+'[ : '+str(len(mesh.edgemax_hystogram[i]))+'/'+str(tot)+' hexes ('+str(len(mesh.edgemax_hystogram[i])/float(tot)*100.)+'%)'+'\n')
+ else:
+ totstat_file.write(str(i)+' ['+str(i*factor+mesh.min_edge_length)+'->'+str((i+1)*factor+mesh.min_edge_length)+'[ : 0/'+str(tot)+' hexes (0%)'+'\n')
+ try:
+ if mesh.dt is not None:
+ totstat_file.write('='*70+'\n')
+ totstat_file.write('STABILITY')
+ totstat_file.write('='*70+'\n')
+ totstat_file.write('time step < '+str(mesh.dt)+'s, for velocity = '+str(mesh.velocity)+'\n')
+ except:
+ pass
+
+ if toclose:
+ totstat_file.close()
+
+ print 'max specfem3d skewness: ',mesh.max_skewness
+ print 'min edge length: ',mesh.min_edge_length
+ return mesh.max_skewness,mesh.min_edge_length
+
+
+def quality_log(tqfile=None):
+ """
+ creation of the quality parameter file
+ """
+ import start as start
+ #
+ #
+ mpiflag,iproc,numproc,mpi = start.start_mpi()
+ #
+ #
+ from hex_metric import SEM_metric_3D
+ #
+ lvol=cubit.parse_cubit_list('volume','all')
+ if len(lvol)!=0:
+ cubit.cmd('quality vol all allmetric ')
+ else:
+ cubit.cmd('quality hex in block all allmetric ')
+ cubit.cmd('list model ')
+ #
+ toclose=True
+ if isinstance(tqfile,file):
+ totstat_file=tqfile
+ elif isinstance(tqfile,str):
+ totstat_file=open(tqfile+'_cubitquality_skewness_proc_'+str(iproc)+'.log','w')
+ else:
+ import sys
+ totstat_file=sys.stdout
+ toclose=False
+
+ mesh=SEM_metric_3D()
+ mesh.check_metric()
+
+ if mesh.max_skewness is not None:
+ mesh.skew_hystogram=mesh.hyst(0,mesh.max_skewness,mesh.skew_hyst)
+ totstat_file.write('-'*70+'\n')
+
+
+
+
+
+ totstat_file.write('='*70+'\n')
+ totstat_file.write('SKEWNESS'+'\n')
+ totstat_file.write('='*70+'\n')
+ if len(mesh.hex_max_skewness) <= 30:
+ totstat_file.write('max = '+str(mesh.max_skewness)+' in hexes '+str(mesh.hex_max_skewness)+'\n')
+ totstat_file.write('(angle -> minimun ='+str(mesh.min_angle)+ ' maximun ='+str(mesh.max_angle)+')'+'\n')
+ else:
+ totstat_file.write('max = '+str(mesh.max_skewness)+' in '+str(len(mesh.hex_max_skewness))+' hexes '+'\n')
+ totstat_file.write('(angle -> minimun ='+str(mesh.min_angle)+' maximun ='+str(mesh.max_angle)+')'+'\n')
+ totstat_file.write('-'*70+'\n')
+ totstat_file.write('skew hystogram')
+ totstat_file.write('-'*70+'\n')
+ tot=0
+ for i in mesh.skew_hystogram.values():
+ tot=tot+len(i)
+ #k=mesh.skew_hystogram.keys()
+ #k.sort()
+ factor=mesh.max_skewness/mesh.nbin
+ for i in range(0,mesh.nbin+1):
+ if mesh.skew_hystogram.has_key(i):
+ if (i+1)*factor <= 1:
+ totstat_file.write(str(i)+' ['+str(i*factor)+'->'+str((i+1)*factor)+'[ : '+str(len(mesh.skew_hystogram[i]))+'/'+str(tot)+' hexes ('+str(len(mesh.skew_hystogram[i])/float(tot)*100.)+'%)'+'\n')
+ else:
+ if (i+1)*factor <= 1:
+ totstat_file.write(str(i)+' ['+str(i*factor)+'->'+str((i+1)*factor)+'[ : 0/'+str(tot)+' hexes (0%)'+'\n')
+ totstat_file.write('-'*70+'\n')
+ ###############################################
+ if mesh.min_edge_length is not None:
+ mesh.edgemin_hystogram=mesh.hyst(mesh.min_edge_length,mesh.max_edge_length,mesh.edgemin_hyst)
+ mesh.edgemax_hystogram=mesh.hyst(mesh.min_edge_length,mesh.max_edge_length,mesh.edgemax_hyst)
+ totstat_file.write('='*70+'\n')
+ totstat_file.write('edge length')
+ totstat_file.write('='*70+'\n')
+ if len(mesh.hex_min_edge_length) <= 30:
+ totstat_file.write('minimum edge length: '+str(mesh.min_edge_length)+ ' in hexes '+str(mesh.hex_min_edge_length)+'\n')
+ else:
+ totstat_file.write('minimum edge length: '+str(mesh.min_edge_length)+ ' in '+str(len(mesh.hex_min_edge_length))+ ' hexes.'+'\n')
+ if len(mesh.hex_max_edge_length) <= 30:
+ totstat_file.write('maximum edge length: '+str(mesh.max_edge_length)+ ' in hexes '+str(mesh.hex_max_edge_length)+'\n')
+ else:
+ totstat_file.write('maximum edge length: '+str(mesh.max_edge_length)+' in '+str(len(mesh.hex_max_edge_length))+ ' hexes.'+'\n')
+ totstat_file.write('-'*70+'\n')
+ totstat_file.write('edge length hystogram')
+ totstat_file.write('-'*70+'\n')
+ factor=(mesh.max_edge_length-mesh.min_edge_length)/mesh.nbin
+ totstat_file.write('minimum edge length'+'\n')
+ tot=0
+ for i in mesh.edgemin_hystogram.values():
+ tot=tot+len(i)
+ #k=mesh.edgemin_hystogram.keys()
+ #k.sort()
+ for i in range(0,mesh.nbin+1):
+ if mesh.edgemin_hystogram.has_key(i):
+ totstat_file.write(str(i)+' ['+str(i*factor+mesh.min_edge_length)+'->'+str((i+1)*factor+mesh.min_edge_length)+'[ : '+str(len(mesh.edgemin_hystogram[i]))+'/'+str(tot)+' hexes ('+str(len(mesh.edgemin_hystogram[i])/float(tot)*100.)+'%)'+'\n')
+ else:
+ totstat_file.write(str(i)+' ['+str(i*factor+mesh.min_edge_length)+'->'+str((i+1)*factor+mesh.min_edge_length)+'[ : 0/'+str(tot)+' hexes (0%)'+'\n')
+ totstat_file.write('-'*70+'\n')
+ totstat_file.write('maximum edge length')
+ tot=0
+ for i in mesh.edgemax_hystogram.values():
+ tot=tot+len(i)
+ #k=mesh.edgemax_hystogram.keys()
+ #k.sort()
+ for i in range(0,mesh.nbin+1):
+ if mesh.edgemax_hystogram.has_key(i):
+ totstat_file.write(str(i)+' ['+str(i*factor+mesh.min_edge_length)+'->'+str((i+1)*factor+mesh.min_edge_length)+'[ : '+str(len(mesh.edgemax_hystogram[i]))+'/'+str(tot)+' hexes ('+str(len(mesh.edgemax_hystogram[i])/float(tot)*100.)+'%)'+'\n')
+ else:
+ totstat_file.write(str(i)+' ['+str(i*factor+mesh.min_edge_length)+'->'+str((i+1)*factor+mesh.min_edge_length)+'[ : 0/'+str(tot)+' hexes (0%)'+'\n')
+ try:
+ if mesh.dt is not None:
+ totstat_file.write('='*70+'\n')
+ totstat_file.write('STABILITY')
+ totstat_file.write('='*70+'\n')
+ totstat_file.write('time step < '+str(mesh.dt)+'s, for velocity = '+str(mesh.velocity)+'\n')
+ except:
+ pass
+
+ if toclose:
+ totstat_file.close()
+
+ print 'max specfem3d skewness: ',mesh.max_skewness
+ print 'min edge length: ',mesh.min_edge_length
+ return mesh.max_skewness,mesh.min_edge_length
\ No newline at end of file
Added: seismo/3D/SPECFEM3D/trunk/CUBIT/read_parameter_cfg.py
===================================================================
--- seismo/3D/SPECFEM3D/trunk/CUBIT/read_parameter_cfg.py (rev 0)
+++ seismo/3D/SPECFEM3D/trunk/CUBIT/read_parameter_cfg.py 2013-09-17 16:14:58 UTC (rev 22796)
@@ -0,0 +1,330 @@
+#############################################################################
+# read_parameter_cfg.py
+# this file is part of GEOCUBIT #
+# #
+# Created by Emanuele Casarotti #
+# Copyright (c) 2008 Istituto Nazionale di Geofisica e Vulcanologia #
+# #
+#############################################################################
+# #
+# GEOCUBIT is free software: you can redistribute it and/or modify #
+# it under the terms of the GNU General Public License as published by #
+# the Free Software Foundation, either version 3 of the License, or #
+# (at your option) any later version. #
+# #
+# GEOCUBIT is distributed in the hope that it will be useful, #
+# but WITHOUT ANY WARRANTY; without even the implied warranty of #
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the #
+# GNU General Public License for more details. #
+# #
+# You should have received a copy of the GNU General Public License #
+# along with GEOCUBIT. If not, see <http://www.gnu.org/licenses/>. #
+# #
+#############################################################################
+import os
+
+def readcfg(filename=None,importmenu=False,mpiflag=False):
+ """
+ read the configuration file, filename is defined in the command line arguments (see menu.py)
+ """
+ if importmenu:
+ import menu as menu
+ cfgname=menu.cfg_name
+ id_proc=menu.id_proc
+ create_plane=menu.create_plane
+ menusurface=menu.surface
+ single=menu.single
+ elif filename:
+ cfgname=filename
+ id_proc=0
+ menu=False
+ create_plane=False
+ menusurface=False
+ single=False
+ else:
+ print 'error: no configuration file'
+ import sys
+ sys.exit()
+ #
+ from utilities import geo2utm #here I can use pyproj but I prefere to include a function in pure python in order to avoid an additional installation
+ #
+ #
+ import ConfigParser
+ config = ConfigParser.ConfigParser()
+ #
+ #
+ def converter(s):
+ if s == 'True':
+ value=True
+ elif s == 'False':
+ value=False
+ elif s == 'None':
+ value= None
+ else:
+ if s.count(',') != 0:
+ value=s.split(',')
+ while value.count(''):
+ value.remove('')
+ else:
+ value=s
+ try:
+ if type(value).__name__ == 'str':
+ if str(value).count('.') != 0:
+ value=float(value)
+ else:
+ value=int(value)
+ else:
+ if str(value).count('.') != 0:
+ value=map(float,value)
+ else:
+ value=map(int,value)
+ except:
+ pass
+ return value
+ #
+ def section_dict(section):
+ dict_o = {}
+ options = config.options(section)
+ for option in options:
+ try:
+ value=converter(config.get(section, option))
+ dict_o[option] = value
+ except:
+ dict_o[option] = None
+ return dict_o
+
+ class attrdict(dict):
+ def __init__(self, *args, **kwargs):
+ dict.__init__(self, *args, **kwargs)
+ self.__dict__ = self
+ def __str__(self):
+ names=[]
+ values=[]
+ for name,value in self.items():
+ names.append(name)
+ values.append(value)
+ 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]
+ arc=o[0][0]
+ return '____'
+
+ #
+ dcfg={}
+ #
+ #CONSTANTS
+ dcfg['osystem'] = 'linux'
+ dcfg['debug_cfg']=False
+ dcfg['checkbound']=False
+ dcfg['top_partitioner'] = 10000
+ dcfg['tres']=0.3 #if n is the vertical component of the normal at a surface pointing horizontally, when -tres < n < tres then the surface is vertical
+ dcfg['precision'] =0.02 #precision for the boundary check (0.02 m)
+ #
+ #INIT
+ dcfg['debug']=True
+ dcfg['cubit_info']="on"
+ dcfg['echo_info']="on"
+ dcfg['jou_info']="on"
+ dcfg['jer_info']="on"
+ dcfg['monitored_cpu']=0
+ dcfg['parallel_import']=True
+ dcfg['save_geometry_cubit']=True
+ dcfg['save_surface_cubit']=False
+ dcfg['save_geometry_paraview'] = False #not implemented
+ dcfg['save_geometry_ACIS'] = False #not implemented
+ dcfg['export_exodus_mesh']=False
+ dcfg['manual_adj'] =False
+ dcfg['play_adj'] =False
+ dcfg['no_adj'] =False
+ dcfg['nx']=False
+ dcfg['ny']=False
+ dcfg['nstep']=False
+ dcfg['localdir_is_globaldir']=True
+ dcfg['refinement_depth']=[]
+ dcfg['scratchdir']=None
+ dcfg['map_meshing_type']='regularmap'
+ dcfg['4sideparallel']=True
+ dcfg["outlinebasin_curve"]=False
+ dcfg["transition_curve"]=False
+ dcfg["faulttrace_curve"]=False
+ dcfg['geological_imprint']=False
+ dcfg['number_processor_xi']=1
+ dcfg['number_processor_eta']=1
+ dcfg['filename']=None
+ dcfg['actual_vertical_interval_top_layer']=1
+ dcfg['coarsening_top_layer']=False
+
+
+
+
+ dcfg['nsurf'] = None
+ if cfgname:
+ config.read(cfgname)
+ sections=['cubit.options','simulation.cpu_parameters','geometry.surfaces','geometry.volumes','geometry.volumes.layercake','geometry.volumes.flatcake','geometry.volumes.partitioner','geometry.partitioner','meshing']
+
+
+ for section in sections:
+ try:
+ d=section_dict(section)
+ dcfg.update(d)
+ except:
+ pass
+
+ if dcfg['nsurf']:
+ surface_name=[]
+ num_x=[]
+ num_y=[]
+ xstep=[]
+ ystep=[]
+ step=[]
+ directionx=[]
+ directiony=[]
+ unit2=[]
+ surf_type=[]
+ delimiter=[]
+ nsurf=int(dcfg['nsurf'])
+ for i in range(1,nsurf+1):
+ section='surface'+str(i)+'.parameters'
+ d=section_dict(section)
+ surface_name.append(d['name'])
+ surf_type.append(d['surf_type'])
+ unit2.append(d['unit_surf'])
+ if d['surf_type'] == 'regular_grid':
+ xstep.append(d['step_along_x'])
+ ystep.append(d['step_along_y'])
+ num_x.append(d['number_point_along_x'])
+ num_y.append(d['number_point_along_y'])
+ elif d['surf_type'] == 'skin':
+ step.append(d['step'])
+ try:
+ delimiter.append(d['delimiter'])
+ except:
+ pass
+ directionx.append(d['directionx'])
+ directiony.append(d['directiony'])
+ dcfg['surface_name']=surface_name
+ dcfg['num_x']=num_x
+ dcfg['num_y']=num_y
+ dcfg['xstep']=xstep
+ dcfg['ystep']=ystep
+ dcfg['step']=step
+ dcfg['directionx']=directionx
+ dcfg['directiony']=directiony
+ dcfg['unit2']=unit2
+ dcfg['surf_type']=surf_type
+ dcfg['delimiter']=delimiter
+
+ try:
+ tres=0
+ xmin,ymin=geo2utm(dcfg['longitude_min'],dcfg['latitude_min'],dcfg['unit'])
+ xmax,ymax=geo2utm(dcfg['longitude_max'],dcfg['latitude_max'],dcfg['unit'])
+ dcfg['xmin']=xmin
+ dcfg['ymin']=ymin
+ dcfg['xmax']=xmax
+ dcfg['ymax']=ymax
+ x1,y1=geo2utm(dcfg['longitude_min'],dcfg['latitude_min'],dcfg['unit'])
+ x2,y2=geo2utm(dcfg['longitude_max'],dcfg['latitude_min'],dcfg['unit'])
+ x3,y3=geo2utm(dcfg['longitude_max'],dcfg['latitude_max'],dcfg['unit'])
+ x4,y4=geo2utm(dcfg['longitude_min'],dcfg['latitude_max'],dcfg['unit'])
+ dcfg['x1_box']=x1
+ dcfg['y1_box']=y1
+ dcfg['x2_box']=x2
+ dcfg['y2_box']=y2
+ dcfg['x3_box']=x3
+ dcfg['y3_box']=y3
+ dcfg['x4_box']=x4
+ dcfg['y4_box']=y4
+ dcfg['tres_boundarydetection']=tres
+ except:
+ pass
+
+ cfg=attrdict(dcfg)
+
+ if menu:
+ try:
+ if cfg.working_dir[-1] == '/': cfg.working_dir=cfg.working_dir[:-1]
+ if cfg.working_dir[0] != '/': cfg.working_dir='./'+cfg.working_dir
+ except:
+ cfg.working_dir=os.getcwd()
+
+ try:
+ if cfg.output_dir[-1] == '/': cfg.output_dir=cfg.output_dir[:-1]
+ if cfg.output_dir[0] != '/': cfg.output_dir='./'+cfg.output_dir
+ except:
+ cfg.output_dir=os.getcwd()
+
+ try:
+ if cfg.SPECFEM3D_output_dir[-1] == '/': cfg.SPECFEM3D_output_dir=cfg.SPECFEM3D_output_dir[:-1]
+ if cfg.SPECFEM3D_output_dir[0] != '/': cfg.SPECFEM3D_output_dir='./'+cfg.SPECFEM3D_output_dir
+ except:
+ cfg.SPECFEM3D_output_dir=os.getcwd()
+
+ cfg.single=single
+
+ if menusurface:
+ cfg.nsurf=1
+ cfg.name=[menu.surface_name]
+ cfg.num_x=[menu.num_x]
+ cfg.num_y=[menu.num_x]
+ cfg.unit=[menu.unit]
+ cfg.surf_type=[menu.surf_type]
+ try:
+ cfg.delimiter=[menu.delimiter]
+ except:
+ cfg.delimiter=[' ']
+ cfg.directionx=[menu.directionx]
+ cfg.directiony=[menu.directiony]
+ else:
+ cfg.SPECFEM3D_output_dir=os.getcwd()
+
+
+ if not cfg.number_processor_eta and cfg.nodes:
+ cfg.number_processor_xi,cfg.number_processor_eta=split(cfg.nodes)
+
+ if isinstance(cfg.filename,str): cfg.filename=[cfg.filename]
+
+ try:
+ cfg.nproc_eta=cfg.number_processor_eta
+ cfg.nproc_xi=cfg.number_processor_xi
+ cfg.cpuy=cfg.number_processor_eta
+ cfg.cpux=cfg.number_processor_xi
+ except:
+ pass
+
+ if create_plane:
+ cfg.x1=map(float,menu.x1.split(','))
+ cfg.x2=map(float,menu.x2.split(','))
+ cfg.x3=map(float,menu.x3.split(','))
+ cfg.x4=map(float,menu.x4.split(','))
+ cfg.unit=menu.unit
+ #
+ if menu:
+ cfg.id_proc=menu.id_proc
+ else:
+ cfg.id_proc=id_proc
+ #
+ try:
+ if isinstance(cfg.tripl,int): cfg.tripl=[cfg.tripl]
+ except:
+ pass
+
+
+ return cfg
+
+
+class getparameter(dict):
+ def __init__(self, *args, **kwargs):
+ dict.__init__(self, *args, **kwargs)
+ self.__dict__ = self
+
+
+def split(x):
+ import math
+ c=int(math.sqrt(x))
+ while math.fmod(x,c):
+ c=c+1
+ return c,x/c
\ No newline at end of file
Modified: seismo/3D/SPECFEM3D/trunk/CUBIT/run_cubit2specfem3d.py
===================================================================
--- seismo/3D/SPECFEM3D/trunk/CUBIT/run_cubit2specfem3d.py 2013-09-17 00:06:13 UTC (rev 22795)
+++ seismo/3D/SPECFEM3D/trunk/CUBIT/run_cubit2specfem3d.py 2013-09-17 16:14:58 UTC (rev 22796)
@@ -1,26 +1,35 @@
-#!python
#!/usr/bin/env python
import cubit
import boundary_definition
import cubit2specfem3d
-
import os
import sys
+# define the name of the CUBIT file to convert to SPECFEM format in the line below
+cubit.cmd('open "large_test_cpml.cub"')
+#cubit.cmd('open "test_cmpl_2layers.cub"')
###### This is boundary_definition.py of GEOCUBIT
+
#..... which extracts the bounding faces and defines them into blocks
-#reload(boundary_definition)
-#boundary_definition.entities=['face']
+
+reload(boundary_definition)
+
+# for CUBIT version 14 and higher (now called TRELIS) you may have to change 'face' to 'Face' in the line below
+boundary_definition.entities=['face']
+
+# this command can run in parallel or not
#boundary_definition.define_bc(boundary_definition.entities,parallel=True)
+boundary_definition.define_bc(boundary_definition.entities,parallel=False)
+#### Export to SPECFEM3D format using cubit2specfem3d.py of GEOCUBIT
-#### Export to SESAME format using cubit2specfem3d.py of GEOCUBIT
os.system('mkdir -p MESH')
reload(cubit2specfem3d)
-cubit2specfem3d.export2SESAME('MESH')
-# all files needed by SCOTCH are now in directory MESH
+cubit2specfem3d.export2SPECFEM3D('MESH')
+# all files needed by SCOTCH are now in directory MESH...
+
Added: seismo/3D/SPECFEM3D/trunk/CUBIT/start.py
===================================================================
--- seismo/3D/SPECFEM3D/trunk/CUBIT/start.py (rev 0)
+++ seismo/3D/SPECFEM3D/trunk/CUBIT/start.py 2013-09-17 16:14:58 UTC (rev 22796)
@@ -0,0 +1,166 @@
+#############################################################################
+# start.py
+# this file is part of GEOCUBIT #
+# #
+# Created by Emanuele Casarotti #
+# Copyright (c) 2008 Istituto Nazionale di Geofisica e Vulcanologia #
+# #
+#############################################################################
+# #
+# GEOCUBIT is free software: you can redistribute it and/or modify #
+# it under the terms of the GNU General Public License as published by #
+# the Free Software Foundation, either version 3 of the License, or #
+# (at your option) any later version. #
+# #
+# GEOCUBIT is distributed in the hope that it will be useful, #
+# but WITHOUT ANY WARRANTY; without even the implied warranty of #
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the #
+# GNU General Public License for more details. #
+# #
+# You should have received a copy of the GNU General Public License #
+# along with GEOCUBIT. If not, see <http://www.gnu.org/licenses/>. #
+# #
+#############################################################################
+#
+#
+#
+#method to call the library
+
+def start_mpi():
+ """
+ start mpi, fakempi mimick the mpi function when the run is serial. The object mpi is based upon pyMPI.
+ it returns: mpiflag,iproc,numproc,mpi
+ where
+ mpiflag is True if parallel mesh is on
+ iproc is the id of the processor (0 if the run is serial)
+ numproc is the number of the processor (1 if the run is serial)
+ mpi is the mpi object
+ """
+ import sys
+ try:
+ import menu as menu
+ iproc=menu.id_proc
+ except:
+ iproc=0
+ try:
+ import mpi
+ numproc=mpi.size
+ mpiflag=True
+ if numproc == 1:
+ mpiflag=False
+ else:
+ iproc=mpi.rank
+ except:
+ class fakempi(object):
+ def __init__(self):
+ self.size=1
+ self.rank=0
+ def barrier(self):
+ pass
+ def allgather(self,value):
+ return value
+ def gather(self,value):
+ return value
+ def bcast(self,value):
+ return value
+ def scatter(self,value):
+ return value
+ def send(self,value,values):
+ return
+ def recv(self,value):
+ return value,value
+ mpi=fakempi()
+ numproc=1
+ mpiflag=False
+ return mpiflag,iproc,numproc,mpi
+
+def start_cubit(init=False):
+ """
+ start cubit, it return the cubit object
+ init argument set the monitotr files
+ """
+ import sys,os
+ try:
+ cubit.silent_cmd('comment')
+ except:
+ try:
+ import cubit
+ cubit.init([""])
+ except:
+ print 'error importing cubit'
+ sys.exit()
+ try:
+ if init:
+ from start import start_cfg,start_mpi
+ cfg=start_cfg()
+ mpiflag,iproc,numproc,mpi = start_mpi()
+ cubit.cmd('set logging on file "'+cfg.working_dir+'/cubit_proc_'+str(iproc)+'.log"')
+ cubit.cmd("set echo off")
+ cubit.cmd("set info off")
+ if iproc == cfg.monitored_cpu:
+ cubit.cmd("record '"+cfg.working_dir+"/monitor_"+str(cfg.monitored_cpu)+".jou'")
+ cubit.cmd("set journal on")
+ cubit.cmd("journal error on")
+ d=cfg.__dict__
+ ks=d.keys()
+ ks.sort()
+ for k in ks:
+ if '__' not in k and '<' not in str(d[k]) and d[k] is not None:
+ txt=str(k)+' -----> '+str(d[k])
+ txt=txt.replace("'","").replace('"','')
+ cubit.cmd('comment "'+txt+'"')
+ else:
+ cubit.cmd("set journal "+cfg.jou_info)
+ cubit.cmd("journal error "+cfg.jer_info)
+ d=cfg.__dict__
+ ks=d.keys()
+ ks.sort()
+ for k in ks:
+ if '__' not in k and '<' not in str(d[k]) and d[k] is not None:
+ txt=str(k)+' -----> '+str(d[k])
+ txt=txt.replace("'","").replace('"','')
+ cubit.cmd('comment "'+txt+'"')
+ cubit.cmd("set echo "+cfg.echo_info)
+ cubit.cmd("set info "+cfg.cubit_info)
+ except:
+ print 'error start cubit'
+ sys.exit()
+ return cubit
+
+def start_cfg(filename=None,importmenu=True):
+ """
+ return the object cfg with the parameters of the mesh
+ """
+ import read_parameter_cfg
+ mpiflag,iproc,numproc,mpi=start_mpi()
+ if filename: importmenu=False
+ cfg=read_parameter_cfg.readcfg(filename=filename,importmenu=importmenu,mpiflag=mpiflag)
+ import os
+ try:
+ os.makedirs(cfg.working_dir)
+ except OSError:
+ pass
+ try:
+ os.makedirs(cfg.output_dir)
+ except OSError:
+ pass
+ try:
+ os.makedirs(cfg.SPECFEM3D_output_dir)
+ except OSError:
+ pass
+
+ return cfg
+
+def start_numpy():
+ """
+ import numpy and check if it is installed
+ """
+ import sys
+ try:
+ import numpy
+ except:
+ print 'error importing numpy, please check if numpy is correctly installed'
+ sys.exit()
+ return numpy
+
+
\ No newline at end of file
Added: seismo/3D/SPECFEM3D/trunk/CUBIT/surfaces.py
===================================================================
--- seismo/3D/SPECFEM3D/trunk/CUBIT/surfaces.py (rev 0)
+++ seismo/3D/SPECFEM3D/trunk/CUBIT/surfaces.py 2013-09-17 16:14:58 UTC (rev 22796)
@@ -0,0 +1,308 @@
+#############################################################################
+# surfaces.py
+# this file is part of GEOCUBIT #
+# #
+# Created by Emanuele Casarotti #
+# Copyright (c) 2008 Istituto Nazionale di Geofisica e Vulcanologia #
+# #
+#############################################################################
+# #
+# GEOCUBIT is free software: you can redistribute it and/or modify #
+# it under the terms of the GNU General Public License as published by #
+# the Free Software Foundation, either version 3 of the License, or #
+# (at your option) any later version. #
+# #
+# GEOCUBIT is distributed in the hope that it will be useful, #
+# but WITHOUT ANY WARRANTY; without even the implied warranty of #
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the #
+# GNU General Public License for more details. #
+# #
+# You should have received a copy of the GNU General Public License #
+# along with GEOCUBIT. If not, see <http://www.gnu.org/licenses/>. #
+# #
+#############################################################################
+try:
+ import start as start
+ cubit = start.start_cubit()
+except:
+ try:
+ import cubit
+ except:
+ print 'error importing cubit, check if cubit is installed'
+ pass
+
+def surfaces(filename=None):
+ """creating the surfaces defined in the parameter files
+ #
+ nsurf = number of surfaces
+ surf_type = list of stype of method for the creation of the surface
+ -- regulare_grid (u and v lines)
+ -- skin (u lines)
+ """
+ #
+ import start as start
+ cfg = start.start_cfg(filename=filename)
+ #
+ #
+ for isurface in range(0,cfg.nsurf):
+ surf_type=cfg.surf_type[isurface]
+ if surf_type == 'regular_grid':
+ surface_regular_grid(isurface,cfgname=filename)
+ elif surf_type == 'skin':
+ surface_skin(isurface,cfgname=filename)
+
+def surface_regular_grid(isurface=0,cfgname=None):
+ """
+ create an acis surface from a regular lon/lat/z grid
+ """
+ import sys,os
+ from math import sqrt
+ from utilities import geo2utm
+ import start as start
+ #
+ #
+ cfg = start.start_cfg(cfgname)
+ numpy = start.start_numpy()
+ #
+ def create_line_u(ind,n,step,data,unit):
+ last_curve_store=cubit.get_last_id("curve")
+ command='create curve spline '
+ for i in range(0,n):
+ if i%step == 0:
+ lon,lat,z=data[i+ind][0],data[i+ind][1],data[i+ind][2]
+ x,y=geo2utm(lon,lat,unit)
+ txt=' Position ' + str(x) +' '+ str(y) +' '+ str(z)
+ command=command+txt
+ #print command
+ cubit.silent_cmd(command)
+ last_curve=cubit.get_last_id("curve")
+ if last_curve != last_curve_store:
+ return last_curve
+ else:
+ return 0
+
+ def create_line_v(ind,n,n2,step,data,unit):
+ last_curve_store=cubit.get_last_id("curve")
+ command='create curve spline '
+ for i in range(0,n):
+ if i%step == 0:
+ lon,lat,z=data[n2*i+ind][0],data[n2*i+ind][1],data[n2*i+ind][2]
+ x,y=geo2utm(lon,lat,unit)
+ txt=' Position ' + str(x) +' '+ str(y) +' '+ str(z)
+ command=command+txt
+ #print command
+ cubit.silent_cmd(command)
+ last_curve=cubit.get_last_id("curve")
+ if last_curve != last_curve_store:
+ return last_curve
+ else:
+ return 0
+ #
+ #
+ cubit.cmd("reset")
+ #
+ position=True
+ #
+ #
+ nu= cfg.num_x[isurface]
+ nv= cfg.num_y[isurface]
+ ustep= cfg.xstep[isurface]
+ vstep= cfg.ystep[isurface]
+ exag=1.
+ unit=cfg.unit2[isurface]
+ #
+ #
+ data=numpy.loadtxt(cfg.surface_name[isurface])
+ if len(data) > 100:
+ command = "set echo off"
+ cubit.cmd(command)
+ command = "set journal off"
+ cubit.cmd(command)
+ #
+ u_curve=[]
+ v_curve=[]
+ #
+ for iv in range(0,nv):
+ if iv%vstep == 0.:
+ u=create_line_u(iv*(nu),nu,ustep,data,unit)
+ u_curve.append(u)
+ for iu in range(0,nu):
+ if iu%ustep == 0.:
+ v=create_line_v(iu,nv,nu,ustep,data,unit)
+ v_curve.append(v)
+ #
+ umax=max(u_curve)
+ umin=min(u_curve)
+ vmax=max(v_curve)
+ vmin=min(v_curve)
+ cubitcommand= 'create surface net u curve '+ str( umin )+' to '+str( umax )+ ' v curve '+ str( vmin )+ ' to '+str( vmax )+' heal'
+ cubit.cmd(cubitcommand)
+ command = "del curve all"
+ cubit.cmd(command)
+ suff=cfg.surface_name[isurface].split('/')
+ command = "save as '"+cfg.working_dir+"/surf_"+suff[-1]+".cub' overwrite"
+ cubit.cmd(command)
+ #
+ #
+ #
+ cubit.cmd("set info "+cfg.cubit_info)
+ cubit.cmd("set echo "+cfg.echo_info)
+ cubit.cmd("set journal "+cfg.jou_info)
+
+def surface_skin(isurface=0,cfgname=None):
+ """
+ create an acis surface interpolating no-intersecting lines
+ """
+ import sys,os
+ from math import sqrt
+ from utilities import geo2utm
+ import start as start
+ #
+ #
+ cubit = start.start_cubit()
+ cfg = start.start_cfg(cfgname)
+ #
+ def define_next_line(directionx,directiony,n,data):
+ ndata=len(data)
+ command=''
+ ind=n
+ try:
+ record=data[ind]
+ except:
+ return False,False
+ try:
+ x,y,z=map(float,record.split())
+ except:
+ return False,False
+ txt=' Position ' + record
+ command=command+txt
+ x_store,y_store,z_store = x,y,z
+ icount=1
+ while True:
+ ind+=1
+ if ind >= ndata: return ind,command
+ record=data[ind]
+ try:
+ x,y,z=map(float,record.split())
+ except:
+ return ind,command
+ dx,dy = x-x_store,y-y_store
+ if directionx == 0 and dy/abs(dy) * directiony >= 0:
+ txt=' Position ' + record
+ command=command+txt
+ icount+=1
+ x_store,y_store,z_store = x,y,z
+ elif directiony == 0 and dx/abs(dx) == directionx :
+ txt=' Position ' + record
+ command=command+txt
+ icount+=1
+ x_store,y_store,z_store = x,y,z
+ else:
+ if icount==1:
+ x,y,z=x_store+1e-4*directionx,y_store+1e-4*directiony,z_store
+ txt=' Position ' +str(x)+ ' '+str(y)+ ' '+str(z)
+ command=command+txt
+ return ind,command
+ def create_line(position):
+ if position:
+ last_curve_store=cubit.get_last_id("curve")
+ command='create curve spline '+position
+ cubit.silent_cmd(command)
+ last_curve=cubit.get_last_id("curve")
+ if last_curve != last_curve_store:
+ return last_curve
+ else:
+ return False
+ else:
+ return False
+
+ command = "reset"
+ cubit.cmd(command)
+ #
+ position=True
+ #
+ try:
+ grdfile = open(cfg.surface_name[isurface], 'r')
+ except:
+ raise NameError, 'No such file or directory: '+ str( cfg.surface_name[isurface] )
+ #
+ directionx=cfg.directionx[isurface]
+ directiony=cfg.directiony[isurface]
+ step=cfg.step[isurface]
+ position=True
+ curveskin=[]
+ count_line=0
+ data=grdfile.read().split('\n')
+ ndata=len(data)
+ n=0
+ #
+ #
+ command = "set echo off"
+ cubit.cmd(command)
+ command = "set journal off"
+ cubit.cmd(command)
+ command = "set info off"
+ cubit.cmd(command)
+ #
+ while position:
+ index,position=define_next_line(directionx,directiony,n,data)
+ if n%step == 0:
+ curve=create_line(position)
+ if curve: curveskin.append(curve)
+ elif n%step != 0 and not position:
+ curve=create_line(position)
+ if curve: curveskin.append(curve)
+ n=index
+ umax=max(curveskin)
+ umin=min(curveskin)
+ print 'create surface skin curve '+ str( umin )+' to '+str( umax )
+ cubitcommand= 'create surface skin curve '+ str( umin )+' to '+str( umax )
+ cubit.cmd(cubitcommand)
+ command = "del curve all"
+ cubit.cmd(command)
+ last_surface=cubit.get_last_id("surface")
+ command = "regularize surf "+str(last_surface)
+ cubit.cmd(command)
+ #
+ suff=cfg.surface_name[isurface].split('/')
+ command = "save as '"+cfg.working_dir+"/surf_"+suff[-1]+".cub' overwrite"
+ cubit.cmd(command)
+ #
+ #
+ #
+ cubit.cmd("set info "+cfg.cubit_info)
+ cubit.cmd("set echo "+cfg.echo_info)
+ cubit.cmd("set journal "+cfg.jou_info)
+
+
+def plane(cfg):
+ import sys,os
+ from math import sqrt
+ from utilities import geo2utm
+ import start as start
+ #
+ #
+ cubit = start.start_cubit()
+ #
+ #
+ command = "reset"
+ cubit.cmd(command)
+ #
+ #
+ for p in [cfg.x1,cfg.x2,cfg.x3,cfg.x4]:
+ x_current,y_current=geo2utm(p[0],p[1],cfg.unit)
+ cubitcommand= 'create vertex '+ str( x_current )+ ' ' + str( y_current) +' '+ str( p[2] )
+ cubit.cmd(cubitcommand)
+ #
+ cubitcommand= 'create surface vertex 1 2 3 4'
+ cubit.cmd(cubitcommand)
+ command = "del vertex all"
+ cubit.cmd(command)
+ command = "save as 'plane.cub' overwrite"
+ cubit.cmd(command)
+ #
+ #
+ #
+ cubit.cmd("set info "+cfg.cubit_info)
+ cubit.cmd("set echo "+cfg.echo_info)
+ cubit.cmd("set journal "+cfg.jou_info)
\ No newline at end of file
Added: seismo/3D/SPECFEM3D/trunk/CUBIT/utilities.py
===================================================================
--- seismo/3D/SPECFEM3D/trunk/CUBIT/utilities.py (rev 0)
+++ seismo/3D/SPECFEM3D/trunk/CUBIT/utilities.py 2013-09-17 16:14:58 UTC (rev 22796)
@@ -0,0 +1,427 @@
+#############################################################################
+# utilities.py
+# this file is part of GEOCUBIT #
+# #
+# Created by Emanuele Casarotti #
+# Copyright (c) 2008 Istituto Nazionale di Geofisica e Vulcanologia #
+# #
+#############################################################################
+# #
+# GEOCUBIT is free software: you can redistribute it and/or modify #
+# it under the terms of the GNU General Public License as published by #
+# the Free Software Foundation, either version 3 of the License, or #
+# (at your option) any later version. #
+# #
+# GEOCUBIT is distributed in the hope that it will be useful, #
+# but WITHOUT ANY WARRANTY; without even the implied warranty of #
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the #
+# GNU General Public License for more details. #
+# #
+# You should have received a copy of the GNU General Public License #
+# along with GEOCUBIT. If not, see <http://www.gnu.org/licenses/>. #
+# #
+#############################################################################
+try:
+ import start as start
+ cubit = start.start_cubit()
+except:
+ try:
+ import cubit
+ except:
+ print 'error importing cubit, check if cubit is installed'
+ pass
+
+
+
+
+def snapshot(name=None,i=0,viewnumber=1):
+ """
+ it takes a snapshot of the figure, following the predefined view position.
+ view 1: vector 1 1 1 z up
+ """
+ if name is None:
+ name='snapshot_'+str(i)
+ i=i+1
+ if viewnumber == 1:
+ command = "at 0"
+ cubit.cmd(command)
+ command = "from 1 1 1"
+ cubit.cmd(command)
+ command = "up 0 0 1"
+ cubit.cmd(command)
+ cubit.cmd('graphics autocenter on')
+ cubit.cmd("zoom reset")
+ command = "hardcopy '"+name+".png' png"
+ cubit.cmd(command)
+ return i
+
+def cubit_error_stop(iproc,command,ner):
+ er=cubit.get_error_count()
+ if er > ner:
+ text='"Proc: '+str(iproc)+' ERROR '+str(command)+' number of error '+str(er)+'/'+str(ner)+'"'
+ cubitcommand = 'comment '+text
+ cubit.cmd(cubitcommand)
+ raise NameError, text
+
+def cubit_error_continue(iproc,command,n_er):
+ er=cubit.get_error_count()
+ if er >= n_er:
+ text='"Proc: '+str(iproc)+' ERROR continue '+str(command)+' number of error '+str(er)+' '+str(n_er)+'"'
+ cubit.cmd('comment '+ text)
+ print 'error: ',text
+ return er
+
+def savemesh(mpiflag,iproc=0,filename=None):
+ import start as start
+ cfg = start.start_cfg(filename=filename)
+ mpiflag,iproc,numproc,mpi = start.start_mpi()
+
+ def runsave(meshfile,iproc,filename=None):
+ import start as start
+ cubit = start.start_cubit()
+ cfg = start.start_cfg(filename=filename)
+ flag=0
+ ner=cubit.get_error_count()
+ cubitcommand= 'save as "'+ cfg.output_dir+'/'+meshfile+'.cub'+ '" overwrite'
+ cubit.cmd(cubitcommand)
+ ner2=cubit.get_error_count()
+ if ner == ner2:
+ cubitcommand= 'export mesh "'+ cfg.output_dir+'/'+meshfile+'.e'+ '" dimension 3 block all overwrite'
+ cubit.cmd(cubitcommand)
+ ner2=cubit.get_error_count()
+ if ner == ner2:
+ flag=1
+ return flag
+
+
+ meshfile='mesh_vol_'+str(iproc)
+
+ flagsaved=0
+ infosave=(iproc,flagsaved)
+
+ mpi.barrier()
+ total_saved=mpi.allgather(flagsaved)
+ if isinstance(total_saved,int): total_saved=[total_saved]
+
+ ind=0
+ saving=True
+ while saving:
+ if len(total_saved) != sum(total_saved):
+ #
+ if not flagsaved:
+ flagsaved=runsave(meshfile,iproc,filename=filename)
+ if flagsaved:
+ infosave=(iproc,flagsaved)
+ if numproc > 1:
+ f=open('mesh_saved'+str(iproc),'w')
+ f.close()
+ mpi.barrier()
+ total_saved=mpi.allgather(flagsaved)
+ if isinstance(total_saved,int): total_saved=[total_saved]
+ ind=ind+1
+ else:
+ saving=False
+ if ind > len(total_saved)+10: saving=False
+ print sum(total_saved),'/',len(total_saved),' saved'
+
+ info_total_saved=mpi.allgather(infosave)
+ if isinstance(info_total_saved,int): info_total_saved=[info_total_saved]
+
+ if iproc==0:
+ f=open('mesh_saving.log','w')
+ f.write('\n'.join(str(x) for x in info_total_saved))
+ f.close()
+
+ f=open(cfg.output_dir+'/'+'blocks_'+str(iproc).zfill(5),'w')
+ blocks=cubit.get_block_id_list()
+
+ for block in blocks:
+ name=cubit.get_exodus_entity_name('block',block)
+ element_count = cubit.get_exodus_element_count(block, "block")
+ nattrib=cubit.get_block_attribute_count(block)
+ attr=[cubit.get_block_attribute_value(block,x) for x in range(0,nattrib)]
+ ty=cubit.get_block_element_type(block)
+ f.write(str(block)+' ; '+name+' ; nattr '+str(nattrib)+' ; '+' '.join(str(x) for x in attr)+' ; '+ty+' '+str(element_count)+'\n')
+ f.close()
+
+ import quality_log
+ f=open(cfg.output_dir+'/'+'quality_'+str(iproc).zfill(5),'w')
+ max_skewness,min_length=quality_log.quality_log(f)
+ f.close()
+
+
+ count_hex=[cubit.get_hex_count()]
+ count_node=[cubit.get_node_count()]
+ max_skew=[(iproc,max_skewness)]
+ min_l=[(iproc,min_length)]
+
+ mpi.barrier()
+ total_min_l=mpi.gather(min_l)
+ total_hex=mpi.gather(count_hex)
+ total_node=mpi.gather(count_node)
+ total_max_skew=mpi.gather(max_skew)
+
+
+ mpi.barrier()
+ if iproc == 0:
+ min_total_min_l=min([ms[1] for ms in total_min_l])
+ max_total_max_skew=max([ms[1] for ms in total_max_skew])
+ sum_total_node=sum(total_node)
+ sum_total_hex=sum(total_hex)
+
+ totstat_file=open(cfg.output_dir+'/totstat.log','w')
+ text='hex total number,node total number,max skew, min length\n'
+ totstat_file.write(text)
+
+ text=str(sum_total_hex)+' , '+str(sum_total_node)+' , '+str(max_total_max_skew)+' , '+str(min_total_min_l)+'\n'
+ totstat_file.write(text)
+
+ totstat_file.write(str(total_max_skew))
+ totstat_file.close()
+
+ print 'meshing process end... proc ',iproc
+
+def importgeometry(geometryfile,iproc=0,filename=None):
+ import start as start
+ cfg = start.start_cfg(filename=filename)
+ mpiflag,iproc,numproc,mpi = start.start_mpi()
+
+ if iproc == 0: print 'importing geometry....'
+ a=['ok from '+str(iproc)]
+
+ mpi.barrier()
+ total_a=mpi.allgather(a)
+ if iproc == 0: print total_a
+
+ def runimport(geometryfile,iproc,filename=None):
+ import start as start
+ cubit = start.start_cubit()
+ cfg = start.start_cfg(filename=filename)
+ file1=cfg.output_dir+'/'+geometryfile
+ cubitcommand= 'open "'+ file1+ '" '
+ cubit.cmd(cubitcommand)
+
+ if cfg.parallel_import:
+ runimport(geometryfile,iproc,filename=filename)
+ else:
+ if iproc == 0:
+ runimport(geometryfile,iproc,filename=filename)
+ for i in range(1,mpi.size):
+ mpi.send('import',i)
+ msg,status=mpi.recv(i)
+ else:
+ msg,status=mpi.recv(0)
+ runimport(geometryfile,iproc,filename=filename)
+ mpi.send('ok'+str(iproc),0)
+
+
+def savesurf(iproc=0):
+ savegeometry(iproc=iproc,surf=True)
+
+###################################################################################### BELOW OK
+
+def load_curves(acis_filename):
+ """
+ load the curves from acis files
+ """
+ import os
+ #
+ #
+ print acis_filename
+ if acis_filename and os.path.exists(acis_filename):
+ tmp_curve=cubit.get_last_id("curve")
+ command = "import acis '"+acis_filename+"'"
+ cubit.cmd(command)
+ tmp_curve_after=cubit.get_last_id("curve")
+ 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'
+ curves=None
+ return [curves]
+
+def project_curves(curves,top_surface):
+ """
+ project curves on surface
+ """
+ if not isinstance(curves,list): curves=curves.split()
+ tmpc=[]
+ for curve in curves:
+ command = "project curve "+str(curve)+" onto surface "+str(top_surface)
+ cubit.cmd(command)
+ tmp_curve_after=cubit.get_last_id("curve")
+ tmpc.append(tmp_curve_after)
+ command = "del curve "+str(curve)
+ cubit.cmd(command)
+ return tmpc
+
+def geo2utm(lon,lat,unit,ellipsoid=23):
+ """conversion geocoodinates from geographical to utm
+
+ usage: x,y=geo2utm(lon,lat,unit,ellipsoid=23)
+
+ dafault ellipsoid is 23 = WGS-84,
+ ellipsoid:
+ 1, "Airy"
+ 2, "Australian National"
+ 3, "Bessel 1841"
+ 4, "Bessel 1841 (Nambia] "
+ 5, "Clarke 1866"
+ 6, "Clarke 1880"
+ 7, "Everest"
+ 8, "Fischer 1960 (Mercury] "
+ 9, "Fischer 1968"
+ 10, "GRS 1967"
+ 11, "GRS 1980"
+ 12, "Helmert 1906"
+ 13, "Hough"
+ 14, "International"
+ 15, "Krassovsky"
+ 16, "Modified Airy"
+ 17, "Modified Everest"
+ 18, "Modified Fischer 1960"
+ 19, "South American 1969"
+ 20, "WGS 60"
+ 21, "WGS 66"
+ 22, "WGS-72"
+ 23, "WGS-84"
+
+ unit: 'geo' if the coordinates of the model (lon,lat) are geographical
+ 'utm' if the coordinates of the model (lon,lat) are utm
+
+ x,y: the function return the easting, northing utm coordinates
+ """
+ import LatLongUTMconversion
+ if unit == 'geo' :
+ (zone, x, y) = LatLongUTMconversion.LLtoUTM(ellipsoid, lat, lon)
+ elif unit == 'utm' :
+ x=lon
+ y=lat
+ return x,y
+
+
+def savegeometry(iproc=0,surf=False,filename=None):
+ import start as start
+ cfg = start.start_cfg(filename=filename)
+ mpiflag,iproc,numproc,mpi = start.start_mpi()
+
+ def runsave(geometryfile,iproc,filename=None):
+ import start as start
+ cubit = start.start_cubit()
+ cfg = start.start_cfg(filename=filename)
+ flag=0
+ ner=cubit.get_error_count()
+ cubitcommand= 'save as "'+ cfg.output_dir+'/'+geometryfile+ '" overwrite'
+ cubit.cmd(cubitcommand)
+ ner2=cubit.get_error_count()
+ if ner == ner2:
+ flag=1
+ return flag
+
+ if surf:
+ geometryfile='surf_vol_'+str(iproc)+'.cub'
+ else:
+ geometryfile='geometry_vol_'+str(iproc)+'.cub'
+
+ flagsaved=0
+ infosave=(iproc,flagsaved)
+
+ mpi.barrier()
+ total_saved=mpi.allgather(flagsaved)
+ if isinstance(total_saved,int): total_saved=[total_saved]
+
+ ind=0
+ saving=True
+ while saving:
+ if len(total_saved) != sum(total_saved):
+ #
+ if not flagsaved:
+ flagsaved=runsave(geometryfile,iproc,filename=filename)
+ if flagsaved:
+ infosave=(iproc,flagsaved)
+ if numproc > 1:
+ f=open('geometry_saved'+str(iproc),'w')
+ f.close()
+ mpi.barrier()
+ total_saved=mpi.allgather(flagsaved)
+ if isinstance(total_saved,int): total_saved=[total_saved]
+ ind=ind+1
+ else:
+ saving=False
+ if ind > len(total_saved)+10: saving=False
+ print sum(total_saved),'/',len(total_saved),' saved'
+
+ info_total_saved=mpi.allgather(infosave)
+ if isinstance(info_total_saved,int): info_total_saved=[info_total_saved]
+
+ if iproc==0:
+ f=open('geometry_saving.log','w')
+ f.write('\n'.join(str(x) for x in info_total_saved))
+ f.close()
+
+
+def get_v_h_list(vol_id_list):
+ """return the lists of the cubit ID of vertical/horizontal surface and vertical/horizontal curves
+ where v/h is defined by the distance of the z normal component from the axis direction
+ the parameter cfg.tres is the threshold as for example if
+ normal[2] >= -tres and normal[2] <= tres
+ then the surface is vertical
+ #
+ usage: surf_or,surf_vertical,list_curve_or,list_curve_vertical,bottom,top = get_v_h_list(list_vol)
+ """
+ #
+ tres=0.3
+
+ try:
+ nvol=len(vol_id_list)
+ except:
+ nvol=1
+ vol_id_list=[vol_id_list]
+ surf_vertical=[]
+ surf_or=[]
+ list_curve_vertical=[]
+ list_curve_or=[]
+ #
+ #
+ for id_vol in vol_id_list:
+ lsurf=cubit.get_relatives("volume",id_vol,"surface")
+ for k in lsurf:
+ normal=cubit.get_surface_normal(k)
+ center_point = cubit.get_center_point("surface", k)
+ if normal[2] >= -1*tres and normal[2] <= tres:
+ surf_vertical.append(k)
+ lcurve=cubit.get_relatives("surface",k,"curve")
+ list_curve_vertical=list_curve_vertical+list(lcurve)
+ else:
+ surf_or.append(k)
+ lcurve=cubit.get_relatives("surface",k,"curve")
+ list_curve_or=list_curve_or+list(lcurve)
+ for x in list_curve_or:
+ try:
+ list_curve_vertical.remove(x)
+ except:
+ pass
+ k=surf_or[0]
+ center_point = cubit.get_center_point("surface", k)[2]
+ center_point_top=center_point
+ center_point_bottom=center_point
+ top=k
+ bottom=k
+ for k in surf_or[1:]:
+ center_point = cubit.get_center_point("surface", k)[2]
+ if center_point > center_point_top:
+ center_point_top=center_point
+ top=k
+ elif center_point < center_point_bottom:
+ center_point_bottom=center_point
+ bottom=k
+ surftop=list(cubit.get_adjacent_surfaces("surface", top))
+ for s in surf_vertical:
+ try:
+ surftop.remove(s)
+ except:
+ pass
+ top=surftop
+ bottom=[bottom]
+ return surf_or,surf_vertical,list_curve_or,list_curve_vertical,bottom,top
+
+
Added: seismo/3D/SPECFEM3D/trunk/CUBIT/volumes.py
===================================================================
--- seismo/3D/SPECFEM3D/trunk/CUBIT/volumes.py (rev 0)
+++ seismo/3D/SPECFEM3D/trunk/CUBIT/volumes.py 2013-09-17 16:14:58 UTC (rev 22796)
@@ -0,0 +1,551 @@
+#############################################################################
+# volumes.py
+# this file is part of GEOCUBIT #
+# #
+# Created by Emanuele Casarotti #
+# Copyright (c) 2008 Istituto Nazionale di Geofisica e Vulcanologia #
+# #
+#############################################################################
+# #
+# GEOCUBIT is free software: you can redistribute it and/or modify #
+# it under the terms of the GNU General Public License as published by #
+# the Free Software Foundation, either version 3 of the License, or #
+# (at your option) any later version. #
+# #
+# GEOCUBIT is distributed in the hope that it will be useful, #
+# but WITHOUT ANY WARRANTY; without even the implied warranty of #
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the #
+# GNU General Public License for more details. #
+# #
+# You should have received a copy of the GNU General Public License #
+# along with GEOCUBIT. If not, see <http://www.gnu.org/licenses/>. #
+# #
+#############################################################################
+try:
+ import start as start
+ cubit = start.start_cubit()
+except:
+ try:
+ import cubit
+ except:
+ print 'error importing cubit, check if cubit is installed'
+ pass
+
+
+def volumes(filename=None):
+ """create the volumes"""
+ import start as start
+ print'volume'
+ cfg = start.start_cfg(filename=filename)
+ #
+ if cfg.volume_type == 'layercake_volume_ascii_regulargrid_regularmap':
+ layercake_volume_ascii_regulargrid_mpiregularmap(filename=filename)
+ elif cfg.volume_type == 'layercake_volume_fromacis_mpiregularmap':
+ layercake_volume_fromacis_mpiregularmap(filename=filename)
+
+def layercake_volume_ascii_regulargrid_mpiregularmap(filename=None):
+ import sys
+ import start as start
+ #
+ mpiflag,iproc,numproc,mpi = start.start_mpi()
+ #
+ numpy = start.start_numpy()
+ cfg = start.start_cfg(filename=filename)
+
+ from utilities import geo2utm, savegeometry,savesurf
+
+ from math import sqrt
+ #
+ try:
+ mpi.barrier()
+ except:
+ pass
+ #
+ #
+ command = "comment '"+"PROC: "+str(iproc)+"/"+str(numproc)+" '"
+ cubit.cmd(command)
+ #fer=open('ERROR_VOLUME_'+str(iproc),'w')
+
+ #
+ if mpiflag:
+ x_slice=numpy.zeros([numproc],int)
+ y_slice=numpy.zeros([numproc],int)
+ for icpuy in range(0,cfg.nproc_eta):
+ for icpux in range (0,cfg.nproc_xi):
+ iprocnum=icpuy*cfg.nproc_xi+icpux
+ #print iprocnum,cfg.nproc_xi,icpux,icpuy,cfg.nproc_xi,cfg.nproc_eta
+ x_slice[iprocnum]=icpux
+ y_slice[iprocnum]=icpuy
+ icpux=x_slice[iproc]
+ icpuy=y_slice[iproc]
+ else:
+ icpuy=int(cfg.id_proc/cfg.nproc_xi)
+ icpux=cfg.id_proc%cfg.nproc_xi
+
+ #
+ if cfg.geometry_format == 'ascii':
+ #for the original surfaces
+ #number of points in the files that describe the topography
+ import local_volume
+ if cfg.localdir_is_globaldir:
+ if iproc == 0 or not mpiflag:
+ coordx_0,coordy_0,elev_0,nx_0,ny_0=local_volume.read_grid(filename)
+ print 'end of reading grd files '+str(nx_0*ny_0)+ ' points'
+ else:
+ pass
+ if iproc == 0 or not mpiflag:
+ coordx=mpi.bcast(coordx_0)
+ else:
+ coordx=mpi.bcast()
+ if iproc == 0 or not mpiflag:
+ coordy=mpi.bcast(coordy_0)
+ else:
+ coordy=mpi.bcast()
+ if iproc == 0 or not mpiflag:
+ elev=mpi.bcast(elev_0)
+ else:
+ elev=mpi.bcast()
+ if iproc == 0 or not mpiflag:
+ nx=mpi.bcast(nx_0)
+ else:
+ nx=mpi.bcast()
+ if iproc == 0 or not mpiflag:
+ ny=mpi.bcast(ny_0)
+ else:
+ ny=mpi.bcast()
+ else:
+ coordx,coordy,elev,nx,ny=local_volume.read_grid(filename)
+ print str(iproc)+ ' end of receving grd files '
+ nx_segment=int(nx/cfg.nproc_xi)+1
+ ny_segment=int(ny/cfg.nproc_eta)+1
+
+ elif cfg.geometry_format=='regmesh': # !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+
+ 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'
+ 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
+ xp=(cfg.xmax-cfg.xmin)/float((nx-1))
+ yp=(cfg.ymax-cfg.ymin)/float((ny-1))
+ #
+ elev=numpy.zeros([nx,ny,cfg.nz],float)
+ coordx=numpy.zeros([nx,ny],float)
+ coordy=numpy.zeros([nx,ny],float)
+ #
+ #
+ xlength=(cfg.xmax-cfg.xmin)/float(cfg.nproc_xi) #length of x slide for chunk
+ ylength=(cfg.ymax-cfg.ymin)/float(cfg.nproc_eta) #length of y slide for chunk
+ nelem_chunk_x=1
+ nelem_chunk_y=1
+ ivxtot=nelem_chunk_x+1
+ 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]
+
+ icoord=0
+ for iy in range(0,ny):
+ for ix in range(0,nx):
+ icoord=icoord+1
+ coordx[ix,iy]=cfg.xmin+xlength*(ix)
+ coordy[ix,iy]=cfg.ymin+ylength*(iy)
+
+ #print coordx,coordy,nx,ny
+ #
+ print 'end of building grid '+str(iproc)
+ print 'number of point: ', len(coordx)
+ #
+ #!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ #for each processor
+ #
+ nxmin_cpu=(nx_segment-1)*(icpux)
+ nymin_cpu=(ny_segment-1)*(icpuy)
+ 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]
+ #
+ #
+ icurve=0
+ isurf=0
+ ivertex=0
+ #
+ #create vertex
+ for inz in range(0,cfg.nz):
+ if cfg.bottomflat and inz == 0: #bottom layer
+ #
+ if cfg.geometry_format == 'ascii':
+ lv=cubit.get_last_id("vertex")
+
+ x_current,y_current=(coordx[nxmin_cpu,nymin_cpu],coordy[nxmin_cpu,nymin_cpu])
+ cubitcommand= 'create vertex '+ str( x_current )+ ' ' + str( y_current) +' '+ str( cfg.depth_bottom )
+ cubit.cmd(cubitcommand)
+ #
+ x_current,y_current=(coordx[nxmin_cpu,nymax_cpu],coordy[nxmin_cpu,nymax_cpu])
+ cubitcommand= 'create vertex '+ str( x_current )+ ' ' + str( y_current) +' '+ str( cfg.depth_bottom )
+ cubit.cmd(cubitcommand)
+ #
+ x_current,y_current=(coordx[nxmax_cpu,nymax_cpu],coordy[nxmax_cpu,nymax_cpu])
+ cubitcommand= 'create vertex '+ str( x_current )+ ' ' + str( y_current) +' '+ str( cfg.depth_bottom )
+ cubit.cmd(cubitcommand)
+ #
+ x_current,y_current=(coordx[nxmax_cpu,nymin_cpu],coordy[nxmax_cpu,nymin_cpu])
+ 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:
+ 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( cfg.depth_bottom )
+ 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( cfg.depth_bottom )
+ 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( 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':
+ 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)
+ #
+ cubitcommand= 'create surface vertex '+str(lv+1)+' '+str(lv+2)+' '+str(lv+3)+' '+str(lv+4)
+ cubit.cmd(cubitcommand)
+ #
+ isurf = isurf + 1
+ elif cfg.geometry_format == 'ascii':
+
+ vertex=[]
+
+ for iy in range(nymin_cpu,nymax_cpu+1):
+ ivx=0
+ for ix in range(nxmin_cpu,nxmax_cpu+1):
+ zvertex=elev[ix,iy,inz]
+ x_current,y_current=(coordx[ix,iy],coordy[ix,iy])
+ #
+ vertex.append(' Position '+ str( x_current ) +' '+ str( y_current )+' '+ str( zvertex ) )
+ #
+ print 'proc',iproc, 'vertex list created....',len(vertex)
+ n=max(nx,ny)
+ uline=[]
+ vline=[]
+ iv=0
+
+ cubit.cmd("set info off")
+ cubit.cmd("set echo off")
+ cubit.cmd("set journal off")
+
+ for iy in range(0,nymax_cpu-nymin_cpu+1):
+ positionx=''
+ for ix in range(0,nxmax_cpu-nxmin_cpu+1):
+ positionx=positionx+vertex[iv]
+ iv=iv+1
+ command='create curve spline '+positionx
+ cubit.cmd(command)
+ #print command
+ uline.append( cubit.get_last_id("curve") )
+ for ix in range(0,nxmax_cpu-nxmin_cpu+1):
+ positiony=''
+ for iy in range(0,nymax_cpu-nymin_cpu+1):
+ positiony=positiony+vertex[ix+iy*(nxmax_cpu-nxmin_cpu+1)]
+ command='create curve spline '+positiony
+ cubit.cmd(command)
+ #print command
+ vline.append( cubit.get_last_id("curve") )
+ #
+ cubit.cmd("set info "+cfg.cubit_info)
+ cubit.cmd("set echo "+cfg.echo_info)
+ cubit.cmd("set journal "+cfg.jou_info)
+ #
+ #
+ print 'proc',iproc, 'lines created....',len(uline),'*',len(vline)
+ umax=max(uline)
+ umin=min(uline)
+ vmax=max(vline)
+ vmin=min(vline)
+ ner=cubit.get_error_count()
+ cubitcommand= 'create surface net u curve '+ str( umin )+' to '+str( umax )+ ' v curve '+ str( vmin )+ ' to '+str( vmax )+' heal'
+ cubit.cmd(cubitcommand)
+ ner2=cubit.get_error_count()
+ if ner == ner2:
+ command = "del curve all"
+ cubit.cmd(command)
+ isurf=isurf+1
+ #
+ else:
+ raise NameError, 'error, check geometry_format, it should be ascii or regmesh' #
+ #
+ cubitcommand= 'del vertex all'
+ cubit.cmd(cubitcommand)
+ if cfg.save_surface_cubit:
+ savegeometry(iproc=iproc,surf=True,filename=filename)
+ #
+ #
+ #!create volume
+ if cfg.osystem == 'macosx':
+ pass
+ elif cfg.osystem == 'linux':
+ for inz in range(1,cfg.nz):
+ ner=cubit.get_error_count()
+ cubitcommand= 'create volume loft surface '+ str( inz+1 )+' '+str( inz )
+ cubit.cmd(cubitcommand)
+ ner2=cubit.get_error_count()
+ isurf=isurf+6
+ if ner == ner2:
+ cubitcommand= 'del surface 1 to '+ str( cfg.nz )
+ cubit.cmd(cubitcommand)
+ list_vol=cubit.parse_cubit_list("volume","all")
+ if len(list_vol) > 1:
+ cubitcommand= 'imprint volume all'
+ cubit.cmd(cubitcommand)
+ cubitcommand= 'merge all'
+ cubit.cmd(cubitcommand)
+ ner=cubit.get_error_count()
+ #cubitcommand= 'composite create curve in vol all'
+ #cubit.cmd(cubitcommand)
+ savegeometry(iproc,filename=filename)
+ if cfg.geological_imprint:
+ curvesname=[cfg.outlinebasin_curve,cfg.transition_curve,cfg.faulttrace_curve]
+ outdir=cfg.working_dir
+ imprint_topography_with_geological_outline(curvesname,outdir)
+ #
+ #
+ cubit.cmd("set info "+cfg.cubit_info)
+ cubit.cmd("set echo "+cfg.echo_info)
+ cubit.cmd("set journal "+cfg.jou_info)
+
+
+def layercake_volume_fromacis_mpiregularmap(filename=None):
+ import sys
+ import start as start
+ #
+ mpiflag,iproc,numproc,mpi = start.start_mpi()
+ #
+ numpy = start.start_numpy()
+ cfg = start.start_cfg(filename=filename)
+ #
+ from utilities import geo2utm, savegeometry
+ #
+ from math import sqrt
+ #
+ try:
+ mpi.barrier()
+ except:
+ pass
+ #
+ #
+ command = "comment '"+"PROC: "+str(iproc)+"/"+str(numproc)+" '"
+ cubit.cmd(command)
+ #
+ #get the limit of the volume considering the cpu
+ def xwebcut(x):
+ command='create planar surface with plane xplane offset '+str(x)
+ cubit.cmd(command)
+ last_surface=cubit.get_last_id("surface")
+ command="webcut volume all tool volume in surf "+str(last_surface)
+ cubit.cmd(command)
+ command="del surf "+str(last_surface)
+ cubit.cmd(command)
+
+ def ywebcut(x):
+ command='create planar surface with plane yplane offset '+str(x)
+ cubit.cmd(command)
+ last_surface=cubit.get_last_id("surface")
+ command="webcut volume all tool volume in surf "+str(last_surface)
+ cubit.cmd(command)
+ command="del surf "+str(last_surface)
+ cubit.cmd(command)
+
+ def translate2zero():
+ ss=cubit.parse_cubit_list('surface','all')
+ box = cubit.get_total_bounding_box("surface", ss)
+ xmin=box[0]
+ ymin=box[3]
+ cubit.cmd('move surface all x '+str(-1*xmin)+' y '+str(-1*ymin))
+ return xmin,ymin
+
+ def translate2original(xmin,ymin):
+ cubit.cmd('move surface all x '+str(xmin)+' y '+str(ymin))
+
+ if mpiflag:
+ x_slice=numpy.zeros([numproc],int)
+ y_slice=numpy.zeros([numproc],int)
+ for icpuy in range(0,cfg.nproc_eta):
+ for icpux in range (0,cfg.nproc_xi):
+ iprocnum=icpuy*cfg.nproc_xi+icpux
+ #print iprocnum,cfg.nproc_xi,icpux,icpuy,cfg.nproc_xi,cfg.nproc_eta
+ x_slice[iprocnum]=icpux
+ y_slice[iprocnum]=icpuy
+ icpux=x_slice[iproc]
+ icpuy=y_slice[iproc]
+ else:
+ icpuy=int(cfg.id_proc/cfg.nproc_xi)
+ icpux=cfg.id_proc%cfg.nproc_xi
+ #
+ ner=cubit.get_error_count()
+ #
+ icurve=0
+ isurf=0
+ ivertex=0
+ #
+ xlength=(cfg.xmax-cfg.xmin)/float(cfg.cpux) #length of x slide for chunk
+ ylength=(cfg.ymax-cfg.ymin)/float(cfg.cpuy) #length of y slide for chunk
+ xmin_cpu=cfg.xmin+(xlength*(icpux))
+ ymin_cpu=cfg.ymin+(ylength*(icpuy))
+ xmax_cpu=xmin_cpu+xlength
+ ymax_cpu=ymin_cpu+ylength
+ #
+ #importing the surfaces
+ for inz in range(cfg.nz-2,-2,-1):
+ if cfg.bottomflat and inz==-1:
+ command = "create planar surface with plane zplane offset "+str(cfg.depth_bottom)
+ cubit.cmd(command)
+ else:
+ command = "import cubit '"+cfg.filename[inz]+"'"
+ cubit.cmd(command)
+
+
+ #translate
+ xmin,ymin=translate2zero()
+ print 'translate ...', -xmin,-ymin
+ xmin_cpu=xmin_cpu-xmin
+ ymin_cpu=ymin_cpu-ymin
+ xmax_cpu=xmax_cpu-xmin
+ ymax_cpu=ymax_cpu-ymin
+
+ ss=cubit.parse_cubit_list('surface','all')
+ box = cubit.get_total_bounding_box("surface", ss)
+ print 'dimension... ', box
+ #cutting the surfaces
+ xwebcut(xmin_cpu)
+ xwebcut(xmax_cpu)
+ ywebcut(ymin_cpu)
+ ywebcut(ymax_cpu)
+ #
+ list_surface_all=cubit.parse_cubit_list("surface","all")
+ #condisidering only the surfaces inside the boundaries
+ dict_surf={}
+ for isurf in list_surface_all:
+ p=cubit.get_center_point("surface",isurf)
+ if p[0] < xmin_cpu or p[0] > xmax_cpu or p[1] > ymax_cpu or p[1] < ymin_cpu:
+ command = "del surf "+str(isurf)
+ cubit.cmd(command)
+ else:
+ dict_surf[str(isurf)]=p[2]
+ z=dict_surf.values()
+ z.sort()
+ list_surf=[]
+ for val in z:
+ isurf=[k for k, v in dict_surf.iteritems() if v == val][0]
+ list_surf.append(int(isurf))
+ #
+
+ #lofting the volume
+ for i,j in zip(list_surf,list_surf[1:]):
+ ner=cubit.get_error_count()
+ cubitcommand= 'create volume loft surface '+ str(i)+' '+str(j)
+ cubit.cmd(cubitcommand)
+ ner2=cubit.get_error_count()
+ #
+ translate2original(xmin,ymin)
+
+
+ if ner == ner2:
+ cubitcommand= 'del surface all'
+ cubit.cmd(cubitcommand)
+ #
+ #
+ #cubitcommand= 'composite create curve in vol all'
+ #cubit.cmd(cubitcommand)
+ list_vol=cubit.parse_cubit_list("volume","all")
+ if len(list_vol) > 1:
+ cubitcommand= 'imprint volume all'
+ cubit.cmd(cubitcommand)
+ #cubit_error_stop(iproc,cubitcommand,ner)
+ #
+ cubitcommand= 'merge all'
+ cubit.cmd(cubitcommand)
+ #
+ savegeometry(iproc,filename=filename)
+ if cfg.geological_imprint:
+ curvesname=[cfg.outlinebasin_curve,cfg.transition_curve,cfg.faulttrace_curve]
+ outdir=cfg.working_dir
+ imprint_topography_with_geological_outline(curvesname,outdir)
+
+
+
+def imprint_topography_with_geological_outline(curvesname,outdir='.'):
+ import sys,os
+ from sets import Set
+ #
+ from utilities import load_curves,project_curves,get_v_h_list
+
+ list_vol=cubit.parse_cubit_list("volume","all")
+ surf_or,surf_vertical,list_curve_or,list_curve_vertical,bottom,top=get_v_h_list(list_vol)
+
+
+
+ outlinebasin_curve=load_curves(curvesname[0])
+ transition_curve=load_curves(curvesname[1])
+ faulttrace_curve=load_curves(curvesname[2])
+
+ curves=[]
+ if outlinebasin_curve: curves=curves+outlinebasin_curve
+ if transition_curve: curves=curves+transition_curve
+ if faulttrace_curve: curves=curves+faulttrace_curve
+
+ if curves:
+ command='imprint tolerant surface '+str(top)+' with curve '+' '.join(str(x) for x in curves)+' merge'
+ cubit.cmd(command)
+
+
+ command = "merge surf all"
+ cubit.cmd(command)
+ command = "compress vol all"
+ cubit.cmd(command)
+ command = "compress surf all"
+ cubit.cmd(command)
+ command = "save as '"+outdirs+"/"+"imprinted_vol_"+str(iproc)+".cub' overwrite"
+ cubit.cmd(command)
+
More information about the CIG-COMMITS
mailing list