[cig-commits] r19254 - in seismo/3D/SPECFEM3D/trunk: examples examples/Mount_StHelens examples/Mount_StHelens/in_data_files examples/homogeneous_halfspace examples/layered_halfspace examples/tomographic_model examples/waterlayered_halfspace src/shared
danielpeter at geodynamics.org
danielpeter at geodynamics.org
Tue Nov 29 22:01:20 PST 2011
Author: danielpeter
Date: 2011-11-29 22:01:19 -0800 (Tue, 29 Nov 2011)
New Revision: 19254
Added:
seismo/3D/SPECFEM3D/trunk/examples/Mount_StHelens/
seismo/3D/SPECFEM3D/trunk/examples/Mount_StHelens/README
seismo/3D/SPECFEM3D/trunk/examples/Mount_StHelens/boundary_definition.py
seismo/3D/SPECFEM3D/trunk/examples/Mount_StHelens/convert_lonlat2utm.pl
seismo/3D/SPECFEM3D/trunk/examples/Mount_StHelens/cubit2specfem3d.py
seismo/3D/SPECFEM3D/trunk/examples/Mount_StHelens/in_data_files/
seismo/3D/SPECFEM3D/trunk/examples/Mount_StHelens/in_data_files/CMTSOLUTION
seismo/3D/SPECFEM3D/trunk/examples/Mount_StHelens/in_data_files/Par_file
seismo/3D/SPECFEM3D/trunk/examples/Mount_StHelens/in_data_files/STATIONS
seismo/3D/SPECFEM3D/trunk/examples/Mount_StHelens/mesh_mount.py
seismo/3D/SPECFEM3D/trunk/examples/Mount_StHelens/process.sh
seismo/3D/SPECFEM3D/trunk/examples/Mount_StHelens/ptopo.mean.utm
seismo/3D/SPECFEM3D/trunk/examples/Mount_StHelens/read_topo.py
seismo/3D/SPECFEM3D/trunk/examples/Mount_StHelens/run_boundary_definition.py
seismo/3D/SPECFEM3D/trunk/examples/Mount_StHelens/run_cubit2specfem3d.py
Modified:
seismo/3D/SPECFEM3D/trunk/examples/homogeneous_halfspace/process.sh
seismo/3D/SPECFEM3D/trunk/examples/layered_halfspace/process.sh
seismo/3D/SPECFEM3D/trunk/examples/tomographic_model/process.sh
seismo/3D/SPECFEM3D/trunk/examples/waterlayered_halfspace/process.sh
seismo/3D/SPECFEM3D/trunk/src/shared/create_movie_shakemap_AVS_DX_GMT.f90
Log:
adds Mount_StHelens CUBIT example; updates some default process.sh scripts
Added: seismo/3D/SPECFEM3D/trunk/examples/Mount_StHelens/README
===================================================================
--- seismo/3D/SPECFEM3D/trunk/examples/Mount_StHelens/README (rev 0)
+++ seismo/3D/SPECFEM3D/trunk/examples/Mount_StHelens/README 2011-11-30 06:01:19 UTC (rev 19254)
@@ -0,0 +1,101 @@
+----------------------------------------------------------------------
+README
+----------------------------------------------------------------------
+
+This example creates a single block model with topography around Mount St.Helens
+using CUBIT, and runs a forward simulation.
+
+step-by-step tutorial:
+
+0. preparing topography data:
+
+ a low resolution 'ptopo.mean.utm' file is provided to run the example,
+ this would be the procedure to create it yourself:
+
+ a) You can get SRTM 90m Digital Elevation Data for a
+ region of interest at: http://srtm.csi.cgiar.org
+
+ For this example, we choose Mount St.Helens as region of interest.
+ Mount St. Helens is located at: 46.197 N 122.186 W
+
+ It is contained in the download: srtm_12_03.zip
+ (see http://srtm.csi.cgiar.org/SELECTION/inputCoord.asp
+ region of interest: Latitude min: 45 N max: 50 N
+ Longitude min: 125 W max: 120 W )
+
+ unzipping the file
+ > unzip srtm_12_03.zip
+ leads to:
+ ..
+ srtm_12_03.tif
+ ..
+
+ b) To convert the tif-file into an lon-lat-elevation format, you can use the
+ package FWTools at: http://fwtools.maptools.org/
+
+ install the package and use their gdal2xyz executable to extract
+ into xyz format:
+ > FWTools-2.0.6/bin_safe/gdal2xyz.py srtm_12_03.tif > srtm_12_03.xyz
+
+ the file 'srtm_12_03.xyz' has now the format: #longitude #latitude #elevation (m)
+ (the file size is ~ 963 MB)
+
+ c) Extract the detailed region of interest:
+
+ To further extract and manipulate the topography data, you can use the
+ package GMT at: http://gmt.soest.hawaii.edu/
+
+ For our purposes, the region of interest will be: -R-122.3/-122.1/46.1/46.3
+ (~23.5 km x 23.5 km )
+
+ Using the blockmean executable, we extract and interpolate the topography data
+ for the detailed region, using an interpolated grid spacing of 0.006 degrees ~ 700 m:
+ > blockmean srtm_12_03.xyz -R-122.3/-122.1/46.1/46.3 -I0.006/0.006 > ptopo.mean.xyz
+
+ c) Since the mesh will need Cartesian coordinates, we convert the topography file
+ from (longitude/latitude/elevation) to UTM (X/Y/Z) coordinates.
+ Mount St.Helens lies in the UTM zone: 10 (T).
+
+ Use the script 'convert_lonlat2utm.pl' provided in this example folder:
+ > ./convert_lonlat2utm.pl ptopo.mean.xyz 10 > ptopo.mean.utm
+ to create a new file 'ptopo.mean.utm' with format: #UTM_X #UTM_Y #Z
+
+
+1. create an STL surface from topography data:
+
+ run script in CUBIT:
+ claro -> Menu "Tools" -> "Play Journal File" ... and select file: "read_topo.py"
+
+ this reads in the input file with name: ptopo.mean.utm
+ and should create a file: topo.stl
+
+ note: CUBIT has limited features when using facets/STL file formats.
+ It is in general preferable to work with ACIS formats.
+
+ However, using an ACIS topography file leads to a slower meshing process.
+ This is why we prefer to use STL file formats here, but rely on
+ a developer command (script was tested using CUBIT version 12.2).
+
+
+2. create mesh:
+
+ - run cubit GUI:
+ > claro
+
+ then run meshing script:
+ claro -> Menu "Tools" -> "Play Journal File" ... and select file: "mesh_mount.py"
+
+ this creates all the mesh files in subdirectory MESH/
+
+
+3. run processing script for a default simulation on 4 CPUs (must have access to mpi):
+
+ - run script 'process.sh' in this example directory:
+ > ./process.sh
+
+ this script handles the following steps:
+ - decomposes mesh files (example creates 4 partitions from mesh files given in MESH/)
+ - runs the database generation
+ - runs a single forward simulation (using the default setup files in in_data_files/)
+
+
Added: seismo/3D/SPECFEM3D/trunk/examples/Mount_StHelens/boundary_definition.py
===================================================================
--- seismo/3D/SPECFEM3D/trunk/examples/Mount_StHelens/boundary_definition.py (rev 0)
+++ seismo/3D/SPECFEM3D/trunk/examples/Mount_StHelens/boundary_definition.py 2011-11-30 06:01:19 UTC (rev 19254)
@@ -0,0 +1,414 @@
+#############################################################################
+# boundary_definition.py #
+# this file is part of GEOCUBIT #
+# #
+# Created by Emanuele Casarotti #
+# Copyright (c) 2008 Istituto Nazionale di Geofisica e Vulcanologia #
+# #
+#############################################################################
+# #
+# GEOCUBIT is free software: you can redistribute it and/or modify #
+# it under the terms of the GNU General Public License as published by #
+# the Free Software Foundation, either version 3 of the License, or #
+# (at your option) any later version. #
+# #
+# GEOCUBIT is distributed in the hope that it will be useful, #
+# but WITHOUT ANY WARRANTY; without even the implied warranty of #
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the #
+# GNU General Public License for more details. #
+# #
+# You should have received a copy of the GNU General Public License #
+# along with GEOCUBIT. If not, see <http://www.gnu.org/licenses/>. #
+# #
+#############################################################################
+
+def define_absorbing_surf():
+ """
+ define the absorbing surfaces for a layered topological box where boundary are surfaces parallel to the axis.
+ it returns absorbing_surf,absorbing_surf_xmin,absorbing_surf_xmax,absorbing_surf_ymin,absorbing_surf_ymax,absorbing_surf_bottom,topo_surf
+ where
+ absorbing_surf is the list of all the absorbing boundary surf
+ absorbing_surf_xmin is the list of the absorbing boundary surfaces that correnspond to x=xmin
+ ...
+ absorbing_surf_bottom is the list of the absorbing boundary surfaces that correspond to z=zmin
+ """
+ try:
+ cubit.cmd('comment')
+ except:
+ try:
+ import cubit
+ cubit.init([""])
+ except:
+ print 'error importing cubit'
+ import sys
+ sys.exit()
+ absorbing_surf=[]
+ absorbing_surf_xmin=[]
+ absorbing_surf_xmax=[]
+ absorbing_surf_ymin=[]
+ absorbing_surf_ymax=[]
+ absorbing_surf_bottom=[]
+ top_surf=[]
+
+
+ list_vol=cubit.parse_cubit_list("volume","all")
+ init_n_vol=len(list_vol)
+ zmax_box=cubit.get_total_bounding_box("volume",list_vol)[7]
+ zmin_box=cubit.get_total_bounding_box("volume",list_vol)[6] #it is the z_min of the box ... box= xmin,xmax,d,ymin,ymax,d,zmin...
+ xmin_box=cubit.get_total_bounding_box("volume",list_vol)[0]
+ xmax_box=cubit.get_total_bounding_box("volume",list_vol)[1]
+ ymin_box=cubit.get_total_bounding_box("volume",list_vol)[3]
+ ymax_box=cubit.get_total_bounding_box("volume",list_vol)[4]
+ list_surf=cubit.parse_cubit_list("surface","all")
+# for k in list_surf:
+# center_point = cubit.get_center_point("surface", k)
+# if abs((center_point[0] - xmin_box)/xmin_box) <= 0.005:
+# absorbing_surf_xmin.append(k)
+# absorbing_surf.append(k)
+# elif abs((center_point[0] - xmax_box)/xmax_box) <= 0.005:
+# absorbing_surf_xmax.append(k)
+# absorbing_surf.append(k)
+# elif abs((center_point[1] - ymin_box)/ymin_box) <= 0.005:
+# absorbing_surf_ymin.append(k)
+# absorbing_surf.append(k)
+# elif abs((center_point[1] - ymax_box)/ymax_box) <= 0.005:
+# absorbing_surf_ymax.append(k)
+# absorbing_surf.append(k)
+# elif abs((center_point[2] - zmin_box)/zmin_box) <= 0.005:
+# absorbing_surf_bottom.append(k)
+# absorbing_surf.append(k)
+# else:
+# sbox=cubit.get_bounding_box('surface',k)
+# dz=abs((sbox[7] - zmax_box)/zmax_box)
+# normal=cubit.get_surface_normal(k)
+# zn=normal[2]
+# dn=abs(zn-1)
+# if dz <= 0.001 and dn < 0.2:
+# top_surf.append(k)
+
+ #box lengths
+ x_len = abs( xmax_box - xmin_box)
+ y_len = abs( ymax_box - ymin_box)
+ z_len = abs( zmax_box - zmin_box)
+
+ print '##boundary box: '
+ print '## x length: ' + str(x_len)
+ print '## y length: ' + str(y_len)
+ print '## z length: ' + str(z_len)
+
+ # tolerance parameters
+ absorbing_surface_distance_tolerance=0.005
+ topographic_surface_distance_tolerance=0.1
+ topographic_surface_normal_tolerance=0.1
+
+ for k in list_surf:
+ center_point = cubit.get_center_point("surface", k)
+ if abs((center_point[0] - xmin_box)/x_len) <= absorbing_surface_distance_tolerance:
+ absorbing_surf_xmin.append(k)
+ absorbing_surf.append(k)
+ elif abs((center_point[0] - xmax_box)/x_len) <= absorbing_surface_distance_tolerance:
+ absorbing_surf_xmax.append(k)
+ absorbing_surf.append(k)
+ elif abs((center_point[1] - ymin_box)/y_len) <= absorbing_surface_distance_tolerance:
+ absorbing_surf_ymin.append(k)
+ absorbing_surf.append(k)
+ elif abs((center_point[1] - ymax_box)/y_len) <= absorbing_surface_distance_tolerance:
+ absorbing_surf_ymax.append(k)
+ absorbing_surf.append(k)
+ elif abs((center_point[2] - zmin_box)/z_len) <= absorbing_surface_distance_tolerance:
+ 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)
+ #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:
+ top_surf.append(k)
+
+ return absorbing_surf,absorbing_surf_xmin,absorbing_surf_xmax,absorbing_surf_ymin,absorbing_surf_ymax,absorbing_surf_bottom,top_surf
+
+def define_absorbing_surf_nopar():
+ """
+ define the absorbing surfaces for a layered topological box where boundary surfaces are not parallel to the axis.
+ it returns absorbing_surf,topo_surf
+ where
+ absorbing_surf is the list of all the absorbing boundary surf
+ """
+ try:
+ cubit.cmd('comment')
+ except:
+ try:
+ import cubit
+ cubit.init([""])
+ except:
+ print 'error importing cubit'
+ import sys
+ sys.exit()
+ from sets import Set
+ def product(*args, **kwds):
+ # product('ABCD', 'xy') --> Ax Ay Bx By Cx Cy Dx Dy
+ # product(range(2), repeat=3) --> 000 001 010 011 100 101 110 111
+ pools = map(tuple, args) * kwds.get('repeat', 1)
+ result = [[]]
+ for pool in pools:
+ result = [x+[y] for x in result for y in pool]
+ return result
+ absorbing_surf=[]
+ absorbing_surf_xmin=[]
+ absorbing_surf_xmax=[]
+ absorbing_surf_ymin=[]
+ absorbing_surf_ymax=[]
+ absorbing_surf_bottom=[]
+ top_surf=[]
+ list_vol=cubit.parse_cubit_list("volume","all")
+ init_n_vol=len(list_vol)
+ zmax_box=cubit.get_total_bounding_box("volume",list_vol)[7]
+ zmin_box=cubit.get_total_bounding_box("volume",list_vol)[6] #it is the z_min of the box ... box= xmin,xmax,d,ymin,ymax,d,zmin...
+ xmin_box=cubit.get_total_bounding_box("volume",list_vol)[0]
+ xmax_box=cubit.get_total_bounding_box("volume",list_vol)[1]
+ ymin_box=cubit.get_total_bounding_box("volume",list_vol)[3]
+ ymax_box=cubit.get_total_bounding_box("volume",list_vol)[4]
+ list_surf=cubit.parse_cubit_list("surface","all")
+ lv=[]
+ for k in list_surf:
+ sbox=cubit.get_bounding_box('surface',k)
+ dzmax=abs((sbox[7] - zmax_box)/zmax_box)
+ dzmin=abs((sbox[6] - zmin_box)/zmin_box)
+ normal=cubit.get_surface_normal(k)
+ zn=normal[2]
+ if dzmax <= 0.001 and zn > 0.7:
+ top_surf.append(k)
+ list_vertex=cubit.get_relatives('surface',k,'vertex')
+ for v in list_vertex:
+ valence=cubit.get_valence(v)
+ if valence <= 4: #valence 3 is a corner, 4 is a vertex between 2 volumes, > 4 is a vertex not in the boundaries
+ lv.append(v)
+ elif dzmin <= 0.001 and zn < -0.7:
+ absorbing_surf.append(k)
+ lp=[]
+ combs=product(lv,lv)
+ for comb in combs:
+ v1=comb[0]
+ v2=comb[1]
+ c=Set(cubit.get_relatives("vertex",v1,"curve")) & Set(cubit.get_relatives("vertex",v2,"curve"))
+ if len(c) == 1:
+ p=cubit.get_center_point("curve",list(c)[0])
+ lp.append(p)
+ for k in list_surf:
+ center_point = cubit.get_center_point("surface", k)
+ for p in lp:
+ if abs((center_point[0] - p[0])/p[0]) <= 0.005 and abs((center_point[1] - p[1])/p[1]) <= 0.005:
+ absorbing_surf.append(k)
+ break
+ return absorbing_surf,top_surf
+
+def define_absorbing_surf_sphere():
+ try:
+ cubit.cmd('comment')
+ except:
+ try:
+ import cubit
+ cubit.init([""])
+ except:
+ print 'error importing cubit'
+ import sys
+ sys.exit()
+ surf=[]
+ list_surf=cubit.parse_cubit_list("surface","all")
+ for s in list_surf:
+ v=cubit.get_relatives('surface',s,'volume')
+ if len(v) == 1:
+ surf.append(s)
+ return surf
+
+def define_block():
+ try:
+ cubit.cmd('comment')
+ except:
+ try:
+ import cubit
+ cubit.init([""])
+ except:
+ print 'error importing cubit'
+ import sys
+ sys.exit()
+ list_vol=cubit.parse_cubit_list("volume","all")
+ init_n_vol=len(list_vol)
+ list_name=map(lambda x: 'vol'+x,map(str,list_vol))
+ return list_vol,list_name
+
+def build_block(vol_list,name):
+ from sets import Set
+ try:
+ cubit.cmd('comment')
+ except:
+ try:
+ import cubit
+ cubit.init([""])
+ except:
+ print 'error importing cubit'
+ import sys
+ sys.exit()
+ block_list=cubit.get_block_id_list()
+ if len(block_list) > 0:
+ id_block=max(block_list)
+ else:
+ id_block=0
+ for v,n in zip(vol_list,name):
+ id_block+=1
+ v_other=Set(vol_list)-Set([v])
+ #command= 'block '+str(id_block)+' hex in node in vol '+str(v)+' except hex in vol '+str(list(v_other))
+ command= 'block '+str(id_block)+' hex in vol '+str(v)
+ command = command.replace("["," ").replace("]"," ")
+ cubit.cmd(command)
+ command = "block "+str(id_block)+" name '"+n+"'"
+ cubit.cmd(command)
+
+def build_block_side(surf_list,name,obj='surface'):
+ try:
+ cubit.cmd('comment')
+ except:
+ try:
+ import cubit
+ cubit.init([""])
+ except:
+ print 'error importing cubit'
+ import sys
+ sys.exit()
+ id_nodeset=cubit.get_next_nodeset_id()
+ id_block=cubit.get_next_block_id()
+
+
+ if obj == 'hex':
+ txt='hex in node in surface'
+ txt1='block '+str(id_block)+ ' '+ txt +' '+str(list(surf_list))
+ txt2="block "+str(id_block)+" name '"+name+"'"
+ txt1=txt1.replace("["," ").replace("]"," ")
+ elif obj == 'node':
+ txt=obj+' in surface'
+ txt1= 'nodeset '+str(id_nodeset)+ ' '+ txt +' '+str(list(surf_list))
+ txt1 = txt1.replace("["," ").replace("]"," ")
+ txt2 = "nodeset "+str(id_nodeset)+" name '"+name+"'"
+ elif obj == 'face' or obj == 'edge':
+ txt=obj+' in surface'
+ txt1= 'block '+str(id_block)+ ' '+ txt +' '+str(list(surf_list))
+ txt1 = txt1.replace("["," ").replace("]"," ")
+ txt2 = "block "+str(id_block)+" name '"+name+"'"
+ else:
+ txt1=''
+ # do not execute: block id might be wrong
+ print "##block "+str(id_block)+" name '"+name+"_notsupported (only hex,face,edge,node)'"
+ txt2=''
+
+
+ cubit.cmd(txt1)
+ cubit.cmd(txt2)
+
+def define_bc(*args,**keys):
+ parallel=keys.get('parallel',True)
+ closed=keys.get('closed',False)
+ if not closed:
+ print "##open region"
+
+ # 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()
+
+ v_list,name_list=define_block()
+ build_block(v_list,name_list)
+ entities=args[0]
+ print entities
+ for entity in entities:
+ print "##entity: "+str(entity)
+
+ # 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:
+ # blocks for each side
+ build_block_side(xmin,entity+'_abs_xmin',obj=entity)
+ build_block_side(xmax,entity+'_abs_xmax',obj=entity)
+ build_block_side(ymin,entity+'_abs_ymin',obj=entity)
+ build_block_side(ymax,entity+'_abs_ymax',obj=entity)
+ build_block_side(bottom,entity+'_abs_bottom',obj=entity)
+
+ # 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)
+
+ else:
+ print "##closed region"
+
+ # 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
+
+
+
+
Added: seismo/3D/SPECFEM3D/trunk/examples/Mount_StHelens/convert_lonlat2utm.pl
===================================================================
--- seismo/3D/SPECFEM3D/trunk/examples/Mount_StHelens/convert_lonlat2utm.pl (rev 0)
+++ seismo/3D/SPECFEM3D/trunk/examples/Mount_StHelens/convert_lonlat2utm.pl 2011-11-30 06:01:19 UTC (rev 19254)
@@ -0,0 +1,285 @@
+#!/usr/bin/perl -w
+
+# extracts longitude/latitude/elevation data from a file
+# and outputs UTM converted X/Y/Z data for a given UTM zone.
+#
+# usage example: ./convert_lonlat2utm.pl ptopo.mean.xyz 10 > ptopo.mean.utm
+
+#use strict;
+use Math::Trig;
+
+sub Usage{
+ print STDERR <<END;
+Usage: convert_lonlat2utm.pl filename zone
+ Requires a file format: lon/lat/elevation
+END
+exit(1)
+}
+
+# gets input arguments
+if (@ARGV != 2) { Usage();}
+$file = $ARGV[0];
+$zone = $ARGV[1];
+
+# checks argument
+if( $zone < 1 or $zone > 60 ) {die("error zone: $zone not UTM zone\n");}
+
+# grab all the locations in file
+open(IN,"$file");
+ at dfile = <IN>;
+close(IN);
+$nlines = @dfile;
+
+for($i=0;$i<$nlines;$i++) {
+ # reads lon/lat coordinates
+ $line = $dfile[$i];
+ ($lon,$lat,$ele) = split(" ",$line);
+
+ # converts to UTM x/y
+ ($x,$y) = geo2utm($lon,$lat,$zone);
+ print "$x \t $y \t $ele \n";
+
+}
+
+#
+#----------------------------------------------------------------------------------------
+#
+
+sub geo2utm {
+ my($lon,$lat,$zone) = @_;
+
+# from utm_geo.f90
+
+#! convert geodetic longitude and latitude to UTM, and back
+#! use iway = ILONGLAT2UTM for long/lat to UTM, IUTM2LONGLAT for UTM to lat/long
+#! a list of UTM zones of the world is available at www.dmap.co.uk/utmworld.htm
+
+# implicit none
+
+# include "constants.h"
+
+#!
+#!-----CAMx v2.03
+#!
+#! UTM_GEO performs UTM to geodetic (long/lat) translation, and back.
+#!
+#! This is a Fortran version of the BASIC program "Transverse Mercator
+#! Conversion", Copyright 1986, Norman J. Berls (Stefan Musarra, 2/94)
+#! Based on algorithm taken from "Map Projections Used by the USGS"
+#! by John P. Snyder, Geological Survey Bulletin 1532, USDI.
+#!
+#! Input/Output arguments:
+#!
+#! rlon Longitude (deg, negative for West)
+#! rlat Latitude (deg)
+#! rx UTM easting (m)
+#! ry UTM northing (m)
+#! UTM_PROJECTION_ZONE UTM zone
+#! iway Conversion type
+#! ILONGLAT2UTM = geodetic to UTM
+#! IUTM2LONGLAT = UTM to geodetic
+#!
+
+ my($UTM_PROJECTION_ZONE,$iway);
+ my($rx,$ry,$rlon,$rlat);
+ my($SUPPRESS_UTM_PROJECTION);
+
+ #my($PI=pi);
+ my($degrad)=pi/180.;
+ my($raddeg)=180.0/pi;
+ my($semimaj)=6378206.4;
+ my($semimin)=6356583.8;
+ my($scfa)=0.9996;
+
+#! some extracts about UTM:
+#!
+#! There are 60 longitudinal projection zones numbered 1 to 60 starting at 180°W.
+#! Each of these zones is 6 degrees wide, apart from a few exceptions around Norway and Svalbard.
+#! There are 20 latitudinal zones spanning the latitudes 80°S to 84°N and denoted
+#! by the letters C to X, ommitting the letter O.
+#! Each of these is 8 degrees south-north, apart from zone X which is 12 degrees south-north.
+#!
+#! To change the UTM zone and the hemisphere in which the
+#! calculations are carried out, need to change the fortran code and recompile. The UTM zone is described
+#! actually by the central meridian of that zone, i.e. the longitude at the midpoint of the zone, 3 degrees
+#! from either zone boundary.
+#! To change hemisphere need to change the "north" variable:
+#! - north=0 for northern hemisphere and
+#! - north=10000000 (10000km) for southern hemisphere. values must be in metres i.e. north=10000000.
+#!
+#! Note that the UTM grids are actually Mercators which
+#! employ the standard UTM scale factor 0.9996 and set the
+#! Easting Origin to 500,000;
+#! the Northing origin in the southern
+#! hemisphere is kept at 0 rather than set to 10,000,000
+#! and this gives a uniform scale across the equator if the
+#! normal convention of selecting the Base Latitude (origin)
+#! at the equator (0 deg.) is followed. Northings are
+#! positive in the northern hemisphere and negative in the
+#! southern hemisphere.
+ my($north)=0.0;
+ my($east)=500000.0;
+
+ my($e2,$e4,$e6,$ep2,$xx,$yy,$dlat,$dlon,$cm,$cmr,$delam);
+ my($f1,$f2,$f3,$f4,$rm,$rn,$t,$c,$a,$e1,$u,$rlat1,$dlat1,$c1,$t1,$rn1,$r1,$d);
+ my($rx_save,$ry_save,$rlon_save,$rlat_save);
+
+ my($IUTM2LONGLAT)=1;
+ my($ILONGLAT2UTM)=2;
+
+#---------------------------------------------------------------
+ $iway = $ILONGLAT2UTM;
+ $UTM_PROJECTION_ZONE = $zone;
+
+ $rlon = $lon;
+ $rlat = $lat;
+#---------------------------------------------------------------
+
+
+#! save original parameters
+ $rlon_save = $rlon;
+ $rlat_save = $rlat;
+ $rx_save = $rx;
+ $ry_save = $ry;
+
+#! define parameters of reference ellipsoid
+ $e2=1.0-($semimin/$semimaj)**2.0;
+ $e4=$e2*$e2;
+ $e6=$e2*$e4;
+ $ep2=$e2/(1.-$e2);
+
+ if ($iway == $IUTM2LONGLAT){
+ $xx = $rx;
+ $yy = $ry;
+ }
+ else{
+ $dlon = $rlon;
+ $dlat = $rlat;
+ }
+
+#!
+#!----- Set Zone parameters
+#!
+ $zone = $UTM_PROJECTION_ZONE;
+ $cm = $zone*6.0 - 183.0; # ! set central meridian for this zone
+ $cmr = $cm*$degrad;
+#!
+#!---- Lat/Lon to UTM conversion
+#!
+ if ($iway == $ILONGLAT2UTM) {
+
+ $rlon = $degrad*$dlon;
+ $rlat = $degrad*$dlat;
+
+ $delam = $dlon - $cm;
+ if ($delam < -180.) {$delam = $delam + 360.;}
+ if ($delam > 180.) {$delam = $delam - 360.;}
+ $delam = $delam*$degrad;
+
+ $f1 = (1. - $e2/4. - 3.*$e4/64. - 5.*$e6/256)*$rlat;
+ $f2 = 3.*$e2/8. + 3.*$e4/32. + 45.*$e6/1024.;
+ $f2 = $f2*sin(2.*$rlat);
+ $f3 = 15.*$e4/256.*45.*$e6/1024.;
+ $f3 = $f3*sin(4.*$rlat);
+ $f4 = 35.*$e6/3072.;
+ $f4 = $f4*sin(6.*$rlat);
+ $rm = $semimaj*($f1 - $f2 + $f3 - $f4);
+ if ($dlat == 90. or $dlat == -90.){
+ $xx = 0.;
+ $yy = $scfa*$rm;
+ }
+ else{
+ $rn = $semimaj/sqrt(1. - $e2*sin($rlat)**2);
+ $t = tan($rlat)**2;
+ $c = $ep2*cos($rlat)**2;
+ $a = cos($rlat)*$delam;
+
+ $f1 = (1. - $t + $c)*$a**3/6.;
+ $f2 = 5. - 18.*$t + $t**2 + 72.*$c - 58.*$ep2;
+ $f2 = $f2*$a**5/120.;
+ $xx = $scfa*$rn*($a + $f1 + $f2);
+ $f1 = $a**2/2.;
+ $f2 = 5. - $t + 9.*$c + 4.*$c**2;
+ $f2 = $f2*$a**4/24.;
+ $f3 = 61. - 58.*$t + $t**2 + 600.*$c - 330.*$ep2;
+ $f3 = $f3*$a**6/720.;
+ $yy = $scfa*($rm + $rn*tan($rlat)*($f1 + $f2 + $f3));
+ }
+ $xx = $xx + $east;
+ $yy = $yy + $north;
+ }
+#!
+#!---- UTM to Lat/Lon conversion
+#!
+ else{
+
+ $xx = $xx - $east;
+ $yy = $yy - $north;
+ $e1 = sqrt(1. - $e2);
+ $e1 = (1. - $e1)/(1. + $e1);
+ $rm = $yy/$scfa;
+ $u = 1. - $e2/4. - 3.*$e4/64. - 5.*$e6/256.;
+ $u = $rm/($semimaj*$u);
+
+ $f1 = 3.*$e1/2. - 27.*$e1**3./32.;
+ $f1 = $f1*sin(2.*$u);
+ $f2 = 21.*$e1**2/16. - 55.*$e1**4/32.;
+ $f2 = $f2*sin(4.*$u);
+ $f3 = 151.*$e1**3./96.;
+ $f3 = $f3*sin(6.*$u);
+ $rlat1 = $u + $f1 + $f2 + $f3;
+ $dlat1 = $rlat1*$raddeg;
+ if ($dlat1 >= 90. or $dlat1 <= -90.){
+ $dlat1 = min($dlat1,90. );
+ $dlat1 = max($dlat1,-90. );
+ $dlon = $cm;
+ }
+ else{
+ $c1 = $ep2*cos($rlat1)**2.;
+ $t1 = tan($rlat1)**2.;
+ $f1 = 1. - $e2*sin($rlat1)**2.;
+ $rn1 = $semimaj/sqrt($f1);
+ $r1 = $semimaj*(1. - $e2)/sqrt($f1**3);
+ $d = $xx/($rn1*$scfa);
+
+ $f1 = $rn1*tan($rlat1)/$r1;
+ $f2 = $d**2/2.;
+ $f3 = 5.*3.*$t1 + 10.*$c1 - 4.*$c1**2 - 9.*$ep2;
+ $f3 = $f3*$d**2*$d**2/24.;
+ $f4 = 61. + 90.*$t1 + 298.*$c1 + 45.*$t1**2. - 252.*$ep2 - 3.*$c1**2;
+ $f4 = $f4*($d**2)**3./720.;
+ $rlat = $rlat1 - $f1*($f2 - $f3 + $f4);
+ $dlat = $rlat*$raddeg;
+
+ $f1 = 1. + 2.*$t1 + $c1;
+ $f1 = $f1*$d**2*$d/6.;
+ $f2 = 5. - 2.*$c1 + 28.*$t1 - 3.*$c1**2 + 8.*$ep2 + 24.*$t1**2.;
+ $f2 = $f2*($d**2)**2*$d/120.;
+ $rlon = $cmr + ($d - $f1 + $f2)/cos($rlat1);
+ $dlon = $rlon*$raddeg;
+ if ($dlon < -180.) {$dlon = $dlon + 360.;}
+ if ($dlon > 180.) {$dlon = $dlon - 360.;}
+ }
+ }
+
+ if ($iway == $IUTM2LONGLAT) {
+ $rlon = $dlon;
+ $rlat = $dlat;
+ $rx = $rx_save;
+ $ry = $ry_save;
+
+ return($rlon,$rlat);
+
+ }
+ else{
+ $rx = $xx;
+ $ry = $yy;
+ $rlon = $rlon_save;
+ $rlat = $rlat_save;
+
+ return($rx,$ry);
+
+ }
+
+
+}
Property changes on: seismo/3D/SPECFEM3D/trunk/examples/Mount_StHelens/convert_lonlat2utm.pl
___________________________________________________________________
Name: svn:executable
+ *
Added: seismo/3D/SPECFEM3D/trunk/examples/Mount_StHelens/cubit2specfem3d.py
===================================================================
--- seismo/3D/SPECFEM3D/trunk/examples/Mount_StHelens/cubit2specfem3d.py (rev 0)
+++ seismo/3D/SPECFEM3D/trunk/examples/Mount_StHelens/cubit2specfem3d.py 2011-11-30 06:01:19 UTC (rev 19254)
@@ -0,0 +1,775 @@
+#!python
+#############################################################################
+# cubit2specfem3d.py #
+# this file is part of GEOCUBIT #
+# #
+# Created by Emanuele Casarotti #
+# Copyright (c) 2008 Istituto Nazionale di Geofisica e Vulcanologia #
+# #
+#############################################################################
+# #
+# GEOCUBIT is free software: you can redistribute it and/or modify #
+# it under the terms of the GNU General Public License as published by #
+# the Free Software Foundation, either version 3 of the License, or #
+# (at your option) any later version. #
+# #
+# GEOCUBIT is distributed in the hope that it will be useful, #
+# but WITHOUT ANY WARRANTY; without even the implied warranty of #
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the #
+# GNU General Public License for more details. #
+# #
+# You should have received a copy of the GNU General Public License #
+# along with GEOCUBIT. If not, see <http://www.gnu.org/licenses/>. #
+# #
+#############################################################################
+#
+#for a complete definition of the format of the mesh in SPECFEM3D_SESAME check the manual (http:/):
+#
+#USAGE
+#
+#############################################################################
+#PREREQUISITE
+#The mesh must be prepared
+# automatically using the module boundary_definition (see boundary_definition.py for more information)
+#or
+# manually following the convention:
+# - each material should have a block defined by:
+# material domain_flag (acoustic/elastic/poroelastic)name,flag of the material (integer),p velocity
+# (or the full description: name, flag, vp, vs, rho, Q ... if not present these last 3 parameters will be
+# interpolated by module mat_parameter)
+# - each mesh should have the block definition for the face on the free_surface (topography),
+# the name of this block must be 'face_topo' or you can change the default name in mesh.topo defined in profile.
+# - each mesh should have the block definition for the faces on the absorbing boundaries,
+# one block for each surface with x=Xmin,x=Xmax,y=Ymin,y=Ymax and z=bottom. The names of
+# the blocks should contain the strings "xmin,xmax,ymin,ymax,bottom"
+#
+#############################################################################
+#RUN
+#In a python script or in the cubit python tab call:
+#
+# export2SESAME(path_exporting_mesh_SPECFEM3D_SESAME)
+#
+#the modele create a python class for the mesh: ex. profile=mesh()
+#and it export the files of the mesh needed by the partitioner of SESAME
+#
+#############################################################################
+#OUTPUT
+#The default output are 11 ASCII files:
+#__________________________________________________________________________________________
+#mesh_name='mesh_file' -> the file that contains the connectity of the all mesh
+# format:
+# number of elements
+# id_elements id_node1 id_node2 id_node3 id_node4 id_node5 id_node6 id_node7 id_node8
+# .....
+#
+#__________________________________________________________________________________________
+##nodecoord_name='nodes_coords_file' -> the file that contains the coordinates of the nodes of the all mesh
+# format:
+# number of nodes
+# id_node x_coordinate y_coordinate z_coordinate
+# .....
+#
+#__________________________________________________________________________________________
+##material_name='materials_file' -> the file that contains the material flag of the elements
+# format:
+# id_element flag
+# .....
+#
+#__________________________________________________________________________________________
+##nummaterial_name='nummaterial_velocity_file' -> table of the material properties
+# format:
+# flag rho vp vs 0 0 #full definition of the properties, flag > 0
+# .....
+# flag 'tomography' file_name #for interpolation with tomography
+# .....
+# flag 'interface' file_name flag_for_the_gll_below_the_interface
+# flag_for_the_gll_above_the_interface #for interpolation with interface
+#__________________________________________________________________________________________
+##absname='absorbing_surface_file' -> this file contains all the face in all the absorbing boundaries
+##absname_local='absorbing_surface_file'+'_xmin' -> this file contains all the face in the
+# absorbing boundary defined by x=Xmin
+##absname_local='absorbing_surface_file'+'_xmax' -> this file contains all the face in the
+# absorbing boundary defined by x=Xmax
+##absname_local='absorbing_surface_file'+'_ymin' -> this file contains all the face in the
+# absorbing boundary defined by y=Ymin
+##absname_local='absorbing_surface_file'+'_ymax' -> this file contains all the face in the
+# absorbing boundary defined by y=Ymax
+##absname_local='absorbing_surface_file'+'_bottom' -> this file contains all the face in the
+# absorbing boundary defined by z=bottom
+# format:
+# number of faces
+# id_(element containg the face) id_node1_face id_node2_face id_node3_face id_node4_face
+# ....
+#
+#__________________________________________________________________________________________
+##freename='free_surface_file' -> file with the hex on the free surface (usually the topography)
+# format:
+# number of faces
+# id_(element containg the face) id_node1_face id_node2_face id_node3_face id_node4_face
+#
+#__________________________________________________________________________________________
+# it is possible save only one (or more) file singularly: for example if you want only the nodecoord_file
+# call the module mesh.nodescoord_write(full path name)
+#
+#############################################################################
+
+import cubit
+
+class mtools(object):
+ """docstring for ciao"""
+ def __init__(self,frequency,list_surf,list_vp):
+ super(mtools, self).__init__()
+ self.frequency = frequency
+ self.list_surf = list_surf
+ self.list_vp = list_vp
+ self.ngll=5
+ self.percent_gll=0.172
+ self.point_wavelength=5
+ def __repr__(self):
+ txt='Meshing for frequency up to '+str(self.frequency)+'Hz\n'
+ for surf,vp in zip(self.list_surf,self.list_vp):
+ txt=txt+'surface '+str(surf)+', vp ='+str(vp)+' -> size '+str(self.freq2meshsize(vp)[0])\
+ +' -> dt '+str(self.freq2meshsize(vp)[0])+'\n'
+ return txt
+ def freq2meshsize(self,vp):
+ velocity=vp*.5
+ self.size=(1/2.5)*velocity/self.frequency*(self.ngll-1)/self.point_wavelength
+ self.dt=.4*self.size/vp*self.percent_gll
+ return self.size,self.dt
+ def mesh_it(self):
+ for surf,vp in zip(self.list_surf,self.list_vp):
+ command = "surface "+str(surf)+" size "+str(self.freq2meshsize(vp)[0])
+ cubit.cmd(command)
+ command = "surface "+str(surf)+ 'scheme pave'
+ cubit.cmd(command)
+ command = "mesh surf "+str(surf)
+ cubit.cmd(command)
+
+class block_tools:
+ def __int__(self):
+ pass
+ def create_blocks(self,mesh_entity,list_entity=None,):
+ if mesh_entity =='surface':
+ txt=' face in surface '
+ elif mesh_entity == 'curve':
+ txt=' edge in curve '
+ elif mesh_entity == 'group':
+ txt=' face in group '
+ if list_entity:
+ if not isinstance(list_entity,list):
+ list_entity=[list_entity]
+ for entity in list_entity:
+ iblock=cubit.get_next_block_id()
+ command = "block "+str(iblock)+ txt +str(entity)
+ cubit.cmd(command)
+ def material_file(self,filename):
+ matfile=open(filename,'w')
+ material=[]
+ for record in matfile:
+ mat_name,vp_str=record.split()
+ vp=float(vp_str)
+ material.append([mat_name,vp])
+ self.material=dict(material)
+ def assign_block_material(self,id_block,mat_name,vp=None):
+ try:
+ material=self.material
+ except:
+ material=None
+ cubit.cmd('block '+str(id_block)+' attribute count 2')
+ cubit.cmd('block '+str(id_block)+' attribute index 1 '+str(id_block))
+ if material:
+ if material.has_key(mat_name):
+ cubit.cmd('block '+str(id_block)+' attribute index 2 '+str(material[mat_name]))
+ print 'block '+str(id_block)+' - material '+mat_name+' - vp '+str(material[mat_name])+' from database'
+ elif vp:
+ cubit.cmd('block '+str(id_block)+' attribute index 2 '+str(vp))
+ print 'block '+str(id_block)+' - material '+mat_name+' - vp '+str(vp)
+ else:
+ print 'assignment impossible: check if '+mat_name+' is in the database or specify vp'
+
+class mesh_tools(block_tools):
+ """Tools for the mesh
+ #########
+ dt,edge_dt,freq,edge_freq=seismic_resolution(edges,velocity,bins_d=None,bins_u=None,sidelist=None,ngll=5,np=8)
+ Given the velocity of a list of edges, seismic_resolution provides the minimum Dt
+ required for the stability condition (and the corrisponding edge).
+ Furthermore, given the number of gll point in the element (ngll) and the number
+ of GLL point for wavelength, it provide the maximum resolved frequency.
+ #########
+ length=edge_length(edge)
+ return the length of a edge
+ #########
+ edge_min,length=edge_min_length(surface)
+ given the cubit id of a surface, it return the edge with minimun length
+ #########
+ """
+ def __int__(self):
+ pass
+ def seismic_resolution(self,edges,velocity,bins_d=None,bins_u=None,sidelist=None):
+ """
+ dt,edge_dt,freq,edge_freq=seismic_resolution(edges,velocity,bins_d=None,bins_u=None,sidelist=None,ngll=5,np=8)
+ Given the velocity of a list of edges, seismic_resolution provides the minimum Dt
+ required for the stability condition (and the corrisponding edge).
+ Furthermore, given the number of gll point in the element (ngll) and the number
+ of GLL point for wavelength, it provide the maximum resolved frequency.
+ """
+ ratiostore=1e10
+ dtstore=1e10
+ edgedtstore=-1
+ edgeratiostore=-1
+ for edge in edges:
+ d=self.edge_length(edge)
+ ratio=(1/2.5)*velocity/d*(self.ngll-1)/self.point_wavelength
+ dt=.4*d/velocity*self.percent_gll
+ if dt<dtstore:
+ dtstore=dt
+ edgedtstore=edge
+ if ratio < ratiostore:
+ ratiostore=ratio
+ edgeratiostore=edge
+ try:
+ for bin_d,bin_u,side in zip(bins_d,bins_u,sidelist):
+ if ratio >= bin_d and ratio < bin_u:
+ command = "sideset "+str(side)+" edge "+str(edge)
+ cubit.cmd(command)
+ #print command
+ break
+ except:
+ pass
+ return dtstore,edgedtstore,ratiostore,edgeratiostore
+ def edge_length(self,edge):
+ """
+ length=edge_length(edge)
+ return the length of a edge
+ """
+ from math import sqrt
+ nodes=cubit.get_connectivity('Edge',edge)
+ x0,y0,z0=cubit.get_nodal_coordinates(nodes[0])
+ x1,y1,z1=cubit.get_nodal_coordinates(nodes[1])
+ d=sqrt((x1-x0)**2+(y1-y0)**2+(z1-z0)**2)
+ return d
+ def edge_min_length(self,surface):
+ """
+ edge_min,length=edge_min_length(surface)
+ given the cubit id of a surface, it return the edge with minimun length
+ """
+ from math import sqrt
+ self.dmin=99999
+ edge_store=0
+ command = "group 'list_edge' add edge in surf "+str(surface)
+ command = command.replace("["," ").replace("]"," ")
+ #print command
+ cubit.cmd(command)
+ group=cubit.get_id_from_name("list_edge")
+ edges=cubit.get_group_edges(group)
+ command = "delete group "+ str(group)
+ cubit.cmd(command)
+ for edge in edges:
+ d=self.edge_length(edge)
+ if d<dmin:
+ self.dmin=d
+ edge_store=edge
+ self.edgemin=edge_store
+ return self.edgemin,self.dmin
+ def normal_check(self,nodes,normal):
+ tres=.2
+ p0=cubit.get_nodal_coordinates(nodes[0])
+ p1=cubit.get_nodal_coordinates(nodes[1])
+ p2=cubit.get_nodal_coordinates(nodes[2])
+ a=[p1[0]-p0[0],p1[1]-p0[1],p1[2]-p0[2]]
+ b=[p2[0]-p1[0],p2[1]-p1[1],p2[2]-p1[2]]
+ axb=[a[1]*b[2] - a[2]*b[1], a[2]*b[0] - a[0]*b[2], a[0]*b[1] - a[1]*b[0]]
+ dot=0.0
+ for i in (0,1,2):
+ dot=dot+axb[i]*normal[i]
+ if dot > 0:
+ return nodes
+ elif dot < 0:
+ return nodes[0],nodes[3],nodes[2],nodes[1]
+ else:
+ print 'error: surface normal, dot=0', axb,normal,dot,p0,p1,p2
+ def mesh_analysis(self,frequency):
+ from sets import Set
+ cubit.cmd('set info off')
+ cubit.cmd('set echo off')
+ cubit.cmd('set journal off')
+ bins_d=[0.0001]+range(0,int(frequency)+1)+[1000]
+ bins_u=bins_d[1:]
+ dt=[]
+ ed_dt=[]
+ r=[]
+ ed_r=[]
+ nstart=cubit.get_next_sideset_id()
+ command = "del sideset all"
+ cubit.cmd(command)
+ for bin_d,bin_u in zip(bins_d,bins_u):
+ nsideset=cubit.get_next_sideset_id()
+ command='create sideset '+str(nsideset)
+ cubit.cmd(command)
+ command = "sideset "+str(nsideset)+ " name "+ "'ratio-["+str(bin_d)+"_"+str(bin_u)+"['"
+ cubit.cmd(command)
+ nend=cubit.get_next_sideset_id()
+ sidelist=range(nstart,nend)
+ for block in self.block_mat:
+ name=cubit.get_exodus_entity_name('block',block)
+ velocity=self.material[name][1]
+ if velocity > 0:
+ faces=cubit.get_block_faces(block)
+ edges=[]
+ for face in faces:
+ es=cubit.get_sub_elements("face", face, 1)
+ edges=edges+list(es)
+ edges=Set(edges)
+ dtstore,edgedtstore,ratiostore,edgeratiostore=self.seismic_resolution(edges,\
+ velocity,bins_d,bins_u,sidelist)
+ dt.append(dtstore)
+ ed_dt.append(edgedtstore)
+ r.append(ratiostore)
+ ed_r.append(edgeratiostore)
+ self.ddt=zip(ed_dt,dt)
+ self.dr=zip(ed_r,r)
+ def sorter(x, y):
+ return cmp(x[1],y[1])
+ self.ddt.sort(sorter)
+ self.dr.sort(sorter)
+ print self.ddt,self.dr
+ print 'Deltat minimum => edge:'+str(self.ddt[0][0])+' dt: '+str(self.ddt[0][1])
+ print 'minimum frequency resolved => edge:'+str(self.dr[0][0])+' frequency: '+str(self.dr[0][1])
+ return self.ddt[0],self.dr[0]
+
+class mesh(object,mesh_tools):
+ def __init__(self):
+ super(mesh, self).__init__()
+ self.mesh_name='mesh_file'
+ self.nodecoord_name='nodes_coords_file'
+ self.material_name='materials_file'
+ self.nummaterial_name='nummaterial_velocity_file'
+ self.absname='absorbing_surface_file'
+ self.freename='free_surface_file'
+ self.recname='STATIONS'
+ self.face='QUAD4'
+ self.face2='SHELL4'
+ self.hex='HEX8'
+ self.edge='BAR2'
+ self.topo='face_topo'
+ self.rec='receivers'
+ self.ngll=5
+ self.percent_gll=0.172
+ self.point_wavelength=5
+ self.block_definition()
+ cubit.cmd('compress')
+ def __repr__(self):
+ pass
+ def block_definition(self):
+ block_flag=[]
+ block_mat=[]
+ block_bc=[]
+ block_bc_flag=[]
+ material={}
+ bc={}
+ blocks=cubit.get_block_id_list()
+ for block in blocks:
+ name=cubit.get_exodus_entity_name('block',block)
+ type=cubit.get_block_element_type(block)
+ print block,name,blocks,type,self.hex,self.face
+ # block has hexahedral elements (HEX8)
+ if type == self.hex:
+ flag=None
+ vel=None
+ vs=None
+ rho=None
+ q=0
+ ani=0
+ # material domain id
+ if name.find("acoustic") >= 0 :
+ imaterial = 1
+ elif name.find("elastic") >= 0 :
+ imaterial = 2
+ elif name.find("poroelastic") >= 0 :
+ imaterial = 3
+ else :
+ imaterial = 0
+ print "block: ",name
+ print " could not find appropriate material for this block..."
+ print ""
+ break
+
+ nattrib=cubit.get_block_attribute_count(block)
+ if nattrib != 0:
+ # material flag:
+ # positive => material properties,
+ # negative => interface/tomography domain
+ flag=int(cubit.get_block_attribute_value(block,0))
+ if flag > 0 and nattrib >= 2:
+ # vp
+ vel=cubit.get_block_attribute_value(block,1)
+ if nattrib >= 3:
+ # vs
+ vs=cubit.get_block_attribute_value(block,2)
+ if nattrib >= 4:
+ #density
+ rho=cubit.get_block_attribute_value(block,3)
+ if nattrib >= 5:
+ #Q_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:
+ #anisotropy_flag
+ ani=cubit.get_block_attribute_value(block,5)
+ elif flag < 0:
+ # velocity model
+ vel=name
+ attrib=cubit.get_block_attribute_value(block,1)
+ if attrib == 1:
+ kind='interface'
+ flag_down=cubit.get_block_attribute_value(block,2)
+ flag_up=cubit.get_block_attribute_value(block,3)
+ elif attrib == 2:
+ kind='tomography'
+ else:
+ flag=block
+ vel,vs,rho,q,ani=(name,0,0,0,0)
+ block_flag.append(int(flag))
+ block_mat.append(block)
+ if flag > 0:
+ par=tuple([imaterial,flag,vel,vs,rho,q,ani])
+ elif flag < 0:
+ if kind=='interface':
+ par=tuple([imaterial,flag,kind,name,flag_down,flag_up])
+ elif kind=='tomography':
+ par=tuple([imaterial,flag,kind,name])
+ elif flag==0:
+ par=tuple([imaterial,flag,name])
+ material[block]=par
+ elif (type == self.face) or (type == self.face2) :
+ # block has surface elements (QUAD4 or SHELL4)
+ block_bc_flag.append(4)
+ block_bc.append(block)
+ bc[block]=4 #face has connectivity = 4
+ if name == self.topo: topography_face=block
+ else:
+ # block elements differ from HEX8/QUAD4/SHELL4
+ print '****************************************'
+ print 'block not properly defined:'
+ print ' name:',name
+ print ' type:',type
+ print
+ print 'please check your block definitions!'
+ print
+ print 'only supported types are:'
+ print ' HEX8 for volumes'
+ print ' QUAD4 for surface'
+ print ' SHELL4 for surface'
+ print '****************************************'
+ continue
+
+ nsets=cubit.get_nodeset_id_list()
+ if len(nsets) == 0: self.receivers=None
+ for nset in nsets:
+ name=cubit.get_exodus_entity_name('nodeset',nset)
+ if name == self.rec:
+ self.receivers=nset
+ else:
+ print 'nodeset '+name+' not defined'
+ self.receivers=None
+ try:
+ self.block_mat=block_mat
+ self.block_flag=block_flag
+ self.block_bc=block_bc
+ self.block_bc_flag=block_bc_flag
+ self.material=material
+ self.bc=bc
+ self.topography=topography_face
+ except:
+ print '****************************************'
+ print 'sorry, no blocks or blocks not properly defined'
+ print block_mat
+ print block_flag
+ print block_bc
+ print block_bc_flag
+ print material
+ print bc
+ print topography
+ print '****************************************'
+ def mat_parameter(self,properties):
+ #note: material property acoustic/elastic/poroelastic are defined by the block's name
+ print "#material properties:"
+ print properties
+ imaterial=properties[0]
+ flag=properties[1]
+ if flag > 0:
+ vel=properties[2]
+ if properties[3] is None and type(vel) != str:
+ # velocity model scales with given vp value
+ if vel >= 30:
+ m2km=1000.
+ else:
+ m2km=1.
+ vp=vel/m2km
+ rho=(1.6612*vp-0.472*vp**2+0.0671*vp**3-0.0043*vp**4+0.000106*vp**4)*m2km
+ txt='%1i %3i %20f %20f %20f %1i %1i\n' % (properties[0],properties[1],rho,vel,vel/(3**.5),0,0)
+ elif type(vel) != str:
+ # velocity model given as vp,vs,rho,..
+ #format nummaterials file: #material_domain_id #material_id #rho #vp #vs #Q_mu #anisotropy_flag
+ txt='%1i %3i %20f %20f %20f %20f %2i\n' % (properties[0],properties[1],properties[4], \
+ properties[2],properties[3],properties[5],properties[6])
+ else:
+ txt='%1i %3i %s \n' % (properties[0],properties[1],properties[2])
+ elif flag < 0:
+ if properties[2] == 'tomography':
+ txt='%1i %3i %s %s\n' % (properties[0],properties[1],properties[2],properties[3])
+ elif properties[2] == 'interface':
+ txt='%1i %3i %s %s %1i %1i\n' % (properties[0],properties[1],properties[2],properties[3],\
+ properties[4],properties[5])
+ return txt
+ def nummaterial_write(self,nummaterial_name):
+ print 'Writing '+nummaterial_name+'.....'
+ nummaterial=open(nummaterial_name,'w')
+ for block in self.block_mat:
+ #name=cubit.get_exodus_entity_name('block',block)
+ nummaterial.write(self.mat_parameter(self.material[block]))
+ nummaterial.close()
+ print 'Ok'
+ def mesh_write(self,mesh_name):
+ meshfile=open(mesh_name,'w')
+ print 'Writing '+mesh_name+'.....'
+ num_elems=cubit.get_hex_count()
+ print ' 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)
+ #print len(hexes)
+ for hexa in hexes:
+ #print hexa
+ nodes=cubit.get_connectivity('Hex',hexa)
+ #nodes=self.jac_check(nodes) #is it valid for 3D? TODO
+ txt=('%10i ')% hexa
+ txt=txt+('%10i %10i %10i %10i %10i %10i %10i %10i\n')% nodes[:]
+ meshfile.write(txt)
+ meshfile.close()
+ print 'Ok'
+ def material_write(self,mat_name):
+ mat=open(mat_name,'w')
+ print 'Writing '+mat_name+'.....'
+ for block,flag in zip(self.block_mat,self.block_flag):
+ hexes=cubit.get_block_hexes(block)
+ for hexa in hexes:
+ mat.write(('%10i %10i\n') % (hexa,flag))
+ mat.close()
+ print 'Ok'
+ def nodescoord_write(self,nodecoord_name):
+ nodecoord=open(nodecoord_name,'w')
+ print 'Writing '+nodecoord_name+'.....'
+ node_list=cubit.parse_cubit_list('node','all')
+ num_nodes=len(node_list)
+ print ' number of nodes:',str(num_nodes)
+ nodecoord.write('%10i\n' % num_nodes)
+ #
+ for node in node_list:
+ x,y,z=cubit.get_nodal_coordinates(node)
+ txt=('%10i %20f %20f %20f\n') % (node,x,y,z)
+ nodecoord.write(txt)
+ nodecoord.close()
+ print 'Ok'
+ def free_write(self,freename=None):
+ # free surface
+ cubit.cmd('set info off')
+ cubit.cmd('set echo off')
+ cubit.cmd('set journal off')
+ from sets import Set
+ normal=(0,0,1)
+ if not freename: freename=self.freename
+ # writes free surface file
+ print 'Writing '+freename+'.....'
+ freehex=open(freename,'w')
+ # searches block definition with name face_topo
+ for block,flag in zip(self.block_bc,self.block_bc_flag):
+ if block == self.topography:
+ name=cubit.get_exodus_entity_name('block',block)
+ print ' block name:',name,'id:',block
+ quads_all=cubit.get_block_faces(block)
+ print ' number of faces = ',len(quads_all)
+ dic_quads_all=dict(zip(quads_all,quads_all))
+ freehex.write('%10i\n' % len(quads_all))
+ list_hex=cubit.parse_cubit_list('hex','all')
+ for h in list_hex:
+ faces=cubit.get_sub_elements('hex',h,2)
+ for f in faces:
+ if dic_quads_all.has_key(f):
+ #print f
+ nodes=cubit.get_connectivity('Face',f)
+ nodes_ok=self.normal_check(nodes,normal)
+ txt='%10i %10i %10i %10i %10i\n' % (h,nodes_ok[0],\
+ nodes_ok[1],nodes_ok[2],nodes_ok[3])
+ freehex.write(txt)
+ freehex.close()
+ print 'Ok'
+ cubit.cmd('set info on')
+ cubit.cmd('set echo on')
+ def abs_write(self,absname=None):
+ # absorbing boundaries
+ import re
+ cubit.cmd('set info off')
+ cubit.cmd('set echo off')
+ cubit.cmd('set journal off')
+ from sets import Set
+ if not absname: absname=self.absname
+ #
+ # loops through all block definitions
+ list_hex=cubit.parse_cubit_list('hex','all')
+ for block,flag in zip(self.block_bc,self.block_bc_flag):
+ if block != self.topography:
+ name=cubit.get_exodus_entity_name('block',block)
+ print name,block
+ absflag=False
+ if re.search('xmin',name):
+ filename=absname+'_xmin'
+ normal=(-1,0,0)
+ elif re.search('xmax',name):
+ filename=absname+'_xmax'
+ normal=(1,0,0)
+ elif re.search('ymin',name):
+ filename=absname+'_ymin'
+ normal=(0,-1,0)
+ elif re.search('ymax',name):
+ filename=absname+'_ymax'
+ normal=(0,1,0)
+ elif re.search('bottom',name):
+ filename=absname+'_bottom'
+ normal=(0,0,-1)
+ elif re.search('abs',name):
+ print " ...face_abs - not used so far..."
+ continue
+ else:
+ continue
+ # opens file
+ print 'Writing '+filename+'.....'
+ abshex_local=open(filename,'w')
+ # gets face elements
+ quads_all=cubit.get_block_faces(block)
+ dic_quads_all=dict(zip(quads_all,quads_all))
+ print ' number of faces = ',len(quads_all)
+ abshex_local.write('%10i\n' % len(quads_all))
+ #command = "group 'list_hex' add hex in face "+str(quads_all)
+ #command = command.replace("["," ").replace("]"," ").replace("("," ").replace(")"," ")
+ #cubit.cmd(command)
+ #group=cubit.get_id_from_name("list_hex")
+ #list_hex=cubit.get_group_hexes(group)
+ #command = "delete group "+ str(group)
+ #cubit.cmd(command)
+ for h in list_hex:
+ faces=cubit.get_sub_elements('hex',h,2)
+ for f in faces:
+ if dic_quads_all.has_key(f):
+ nodes=cubit.get_connectivity('Face',f)
+ if not absflag:
+ # checks with specified normal
+ nodes_ok=self.normal_check(nodes,normal)
+ txt='%10i %10i %10i %10i %10i\n' % (h,nodes_ok[0],\
+ nodes_ok[1],nodes_ok[2],nodes_ok[3])
+ else:
+ txt='%10i %10i %10i %10i %10i\n' % (h,nodes[0],\
+ nodes[1],nodes[2],nodes[3])
+ abshex_local.write(txt)
+ # closes file
+ abshex_local.close()
+ print 'Ok'
+ cubit.cmd('set info on')
+ cubit.cmd('set echo on')
+ def surface_write(self,pathdir=None):
+ # optional surfaces, e.g. moho_surface
+ # should be created like e.g.:
+ # > block 10 face in surface 2
+ # > block 10 name 'moho_surface'
+ import re
+ from sets import Set
+ for block in self.block_bc :
+ if block != self.topography:
+ name=cubit.get_exodus_entity_name('block',block)
+ # skips block names like face_abs**, face_topo**
+ if re.search('abs',name):
+ continue
+ elif re.search('topo',name):
+ continue
+ elif re.search('surface',name):
+ filename=pathdir+name+'_file'
+ else:
+ continue
+ # gets face elements
+ print ' surface block name: ',name,'id: ',block
+ quads_all=cubit.get_block_faces(block)
+ print ' face = ',len(quads_all)
+ if len(quads_all) == 0 :
+ continue
+ # writes out surface infos to file
+ print 'Writing '+filename+'.....'
+ surfhex_local=open(filename,'w')
+ dic_quads_all=dict(zip(quads_all,quads_all))
+ # writes number of surface elements
+ surfhex_local.write('%10i\n' % len(quads_all))
+ # writes out element node ids
+ list_hex=cubit.parse_cubit_list('hex','all')
+ for h in list_hex:
+ faces=cubit.get_sub_elements('hex',h,2)
+ for f in faces:
+ if dic_quads_all.has_key(f):
+ nodes=cubit.get_connectivity('Face',f)
+ txt='%10i %10i %10i %10i %10i\n' % (h,nodes[0],\
+ nodes[1],nodes[2],nodes[3])
+ surfhex_local.write(txt)
+ # closes file
+ surfhex_local.close()
+ print 'Ok'
+ def rec_write(self,recname):
+ print 'Writing '+self.recname+'.....'
+ recfile=open(self.recname,'w')
+ nodes=cubit.get_nodeset_nodes(self.receivers)
+ for i,n in enumerate(nodes):
+ x,y,z=cubit.get_nodal_coordinates(n)
+ recfile.write('ST%i XX %20f %20f 0.0 0.0 \n' % (i,x,z))
+ recfile.close()
+ print 'Ok'
+ def write(self,path=''):
+ cubit.cmd('set info off')
+ cubit.cmd('set echo off')
+ cubit.cmd('set journal off')
+ if len(path) != 0:
+ if path[-1] != '/': path=path+'/'
+ # mesh file
+ self.mesh_write(path+self.mesh_name)
+ # mesh material
+ self.material_write(path+self.material_name)
+ # mesh coordinates
+ self.nodescoord_write(path+self.nodecoord_name)
+ # material definitions
+ self.nummaterial_write(path+self.nummaterial_name)
+ # free surface: face_top
+ self.free_write(path+self.freename)
+ # absorbing surfaces: abs_***
+ self.abs_write(path+self.absname)
+ # any other surfaces: ***surface***
+ self.surface_write(path)
+ # receivers
+ if self.receivers: self.rec_write(path+self.recname)
+ cubit.cmd('set info on')
+ cubit.cmd('set echo on')
+
+def export2SESAME(path_exporting_mesh_SPECFEM3D_SESAME):
+ cubit.cmd('set info on')
+ cubit.cmd('set echo on')
+ sem_mesh=mesh()
+ sem_mesh.write(path=path_exporting_mesh_SPECFEM3D_SESAME)
+
+
+if __name__ == '__main__':
+ path='MESH/'
+ export2SESAME(path)
+
+# call by:
+# import cubit2specfem3d
+# cubit2specfem3d.export2SESAME('MESH')
Added: seismo/3D/SPECFEM3D/trunk/examples/Mount_StHelens/in_data_files/CMTSOLUTION
===================================================================
--- seismo/3D/SPECFEM3D/trunk/examples/Mount_StHelens/in_data_files/CMTSOLUTION (rev 0)
+++ seismo/3D/SPECFEM3D/trunk/examples/Mount_StHelens/in_data_files/CMTSOLUTION 2011-11-30 06:01:19 UTC (rev 19254)
@@ -0,0 +1,13 @@
+PDE 1999 01 01 00 00 00.00 67000 67000 -25000 4.2 4.2 hom_explosion
+event name: hom_explosion
+time shift: 0.0000
+half duration: 0.5
+latitude: 46.197
+longitude: -122.186
+depth: 5.0
+Mrr: 1.000000e+23
+Mtt: 1.000000e+23
+Mpp: 1.000000e+23
+Mrt: 0.000000
+Mrp: 0.000000
+Mtp: 0.000000
Added: seismo/3D/SPECFEM3D/trunk/examples/Mount_StHelens/in_data_files/Par_file
===================================================================
--- seismo/3D/SPECFEM3D/trunk/examples/Mount_StHelens/in_data_files/Par_file (rev 0)
+++ seismo/3D/SPECFEM3D/trunk/examples/Mount_StHelens/in_data_files/Par_file 2011-11-30 06:01:19 UTC (rev 19254)
@@ -0,0 +1,53 @@
+
+# forward or adjoint simulation
+SIMULATION_TYPE = 1 # 1 = forward, 2 = adjoint, 3 = both simultaneously
+NOISE_TOMOGRAPHY = 0 # 0 = earthquake simulation, 1/2/3 = three steps in noise simulation
+SAVE_FORWARD = .false.
+
+# UTM projection parameters
+UTM_PROJECTION_ZONE = 10
+SUPPRESS_UTM_PROJECTION = .false.
+
+# number of MPI processors
+NPROC = 4
+
+# time step parameters
+NSTEP = 2500
+DT = 0.005d0
+
+# parameters describing the model
+OCEANS = .false.
+TOPOGRAPHY = .false.
+ATTENUATION = .false.
+USE_OLSEN_ATTENUATION = .false.
+ANISOTROPY = .false.
+
+# absorbing boundary conditions for a regional simulation
+ABSORBING_CONDITIONS = .true.
+
+# save AVS or OpenDX movies
+MOVIE_SURFACE = .false.
+MOVIE_VOLUME = .false.
+NTSTEP_BETWEEN_FRAMES = 200
+CREATE_SHAKEMAP = .false.
+SAVE_DISPLACEMENT = .false.
+USE_HIGHRES_FOR_MOVIES = .false.
+HDUR_MOVIE = 0.0
+
+# save AVS or OpenDX mesh files to check the mesh
+SAVE_MESH_FILES = .true.
+
+# path to store the local database file on each node
+LOCAL_PATH = ../in_out_files/DATABASES_MPI
+
+# interval at which we output time step info and max of norm of displacement
+NTSTEP_BETWEEN_OUTPUT_INFO = 500
+
+# interval in time steps for writing of seismograms
+NTSTEP_BETWEEN_OUTPUT_SEISMOS = 10000
+
+# interval in time steps for reading adjoint traces
+NTSTEP_BETWEEN_READ_ADJSRC = 0 # 0 = read the whole adjoint sources at the same time
+
+# print source time function
+PRINT_SOURCE_TIME_FUNCTION = .false.
Added: seismo/3D/SPECFEM3D/trunk/examples/Mount_StHelens/in_data_files/STATIONS
===================================================================
--- seismo/3D/SPECFEM3D/trunk/examples/Mount_StHelens/in_data_files/STATIONS (rev 0)
+++ seismo/3D/SPECFEM3D/trunk/examples/Mount_StHelens/in_data_files/STATIONS 2011-11-30 06:01:19 UTC (rev 19254)
@@ -0,0 +1,8 @@
+X10 DB 46.197 -122.186 0.0 0.0
+X20 DB 46.2 -122.2 0.0 0.0
+X30 DB 46.29 -122.2 0.0 0.0
+X40 DB 46.11 -122.2 0.0 0.0
+Y10 DB 46.197 -122.186 0.0 0.0
+Y20 DB 46.2 -122.2 0.0 0.0
+Y30 DB 46.2 -122.29 0.0 0.0
+Y40 DB 46.2 -122.11 0.0 0.0
Added: seismo/3D/SPECFEM3D/trunk/examples/Mount_StHelens/mesh_mount.py
===================================================================
--- seismo/3D/SPECFEM3D/trunk/examples/Mount_StHelens/mesh_mount.py (rev 0)
+++ seismo/3D/SPECFEM3D/trunk/examples/Mount_StHelens/mesh_mount.py 2011-11-30 06:01:19 UTC (rev 19254)
@@ -0,0 +1,162 @@
+#!/usr/bin/env python
+
+import cubit
+import boundary_definition
+import cubit2specfem3d
+
+import os
+import sys
+import os.path
+
+cubit.cmd('reset')
+os.system('rm -f topo_brick.stl topo_vol.stl topo_vol2.stl')
+
+print "running meshing script..."
+print ""
+print "note: this script uses topography surface in STL format"
+print " meshing will take around 2 min"
+print ""
+
+# note: this is a workaround to use STL file formats rather than ACIS formats.
+# for our purpose to create a simple mesh, this STL formats are faster
+
+#############################################################
+#
+# 1. step: creates temporary brick volume in STL format
+#
+#############################################################
+
+# creates temporary brick volume and exports in STL format
+# new brick volume (depth will become 1/2 * 20,000 = 10,000 m)
+cubit.cmd('brick x 15000 y 22000 z 20000')
+# moves volume to UTM coordinates of topography surface
+cubit.cmd('volume 1 move x 561738. y 5116370. z 0 ')
+# saves as STL body
+cubit.cmd('export stl ascii "topo_brick.stl" overwrite')
+# clear
+cubit.cmd('reset')
+#checks if new file available
+if not os.path.exists("topo_brick.stl"):
+ print ""
+ print "error creating new STL file topo_brick.stl, please check manually..."
+ print ""
+ cubit.cmd('pause')
+
+#############################################################
+#
+# 2. step: imports topography surface and brick volume in STL format
+#
+#############################################################
+
+# topography surface
+if os.path.exists("topo.stl"):
+ print "opening existing topography surface"
+ # previously run, just reopen the cubit file
+ cubit.cmd('import stl "topo.stl" merge stitch')
+else:
+ # topo surface doesn't exist yet, this creates it:
+ print "reading in topography surface"
+ # reads in topography points and creates sheet surface
+ execfile("./read_topo.py")
+ # clear
+ cubit.cmd('reset')
+ # now reopen the cubit file
+ cubit.cmd('import stl "topo.stl" merge stitch')
+# re-opens brick in STL format
+cubit.cmd('import stl "topo_brick.stl" merge stitch')
+# imprints topography surfaces into brick volume, creates connected surfaces
+# note: this is a develop feature
+cubit.cmd('set developer commands on')
+cubit.cmd('imprint all')
+# exports as STL only surfaces which will create new volume
+cubit.cmd('export stl ascii "topo_vol.stl" surface 9 5 6 8 11 13 overwrite')
+# clear
+cubit.cmd('reset')
+#checks if new file available
+if not os.path.exists("topo_vol.stl"):
+ print ""
+ print "error creating new STL file topo_vol.stl, please check manually..."
+ print ""
+ cubit.cmd('pause')
+
+#############################################################
+#
+# 3. step: manipulate STL file to create a single volume
+#
+#############################################################
+
+os.system('awk \'BEGIN{print \"solid Body_1\";}{if($0 !~ /solid/) print $0;}END{print \"endsolid Body_1\";}\' topo_vol.stl > topo_vol2.stl')
+#checks if new file available
+if not os.path.exists("topo_vol2.stl"):
+ print ""
+ print "error creating new STL file topo_vol2.stl, please check manually..."
+ print ""
+ cubit.cmd('pause')
+
+#############################################################
+#
+# 4. step: import STL volume and create mesh
+#
+#############################################################
+
+# re-opens new volume in STL format
+cubit.cmd('import stl "topo_vol2.stl" merge stitch')
+# Meshing the volumes
+elementsize = 2000.0
+cubit.cmd('volume all size '+str(elementsize))
+# sets meshing type
+# uses a sweep algorithm in vertical (Z) direction
+cubit.cmd('volume all scheme sweep Vector 0 0 1')
+# initial coarse mesh
+cubit.cmd('mesh volume all')
+# draw/update mesh lines for visualization
+# this will draw also the tripling layer mesh lines
+cubit.cmd('draw volume all')
+# optional smoothing to improve mesh quality (takes up to 5 min)
+cubit.cmd('volume all smooth scheme condition number beta 1.0 cpu 5')
+cubit.cmd('smooth volume all')
+# refines global mesh
+cubit.cmd('refine volume 1 numsplit 1')
+cubit.cmd('draw volume all')
+# optional smoothing
+cubit.cmd('smooth volume all')
+# refines elements at topography surface
+cubit.cmd('refine surface 2 numsplit 1 bias 1.0 depth 1')
+cubit.cmd('draw volume all')
+# optional smoothing to improve mesh quality (takes up all 10 min if beta is 2.0)
+cubit.cmd('volume all smooth scheme condition number beta 2.3 cpu 10')
+cubit.cmd('smooth volume all')
+cubit.cmd('draw volume all')
+
+#### End of meshing
+
+###### This is boundary_definition.py of GEOCUBIT
+#..... which extracts the bounding faces and defines them into blocks
+boundary_definition.entities=['face']
+boundary_definition.define_bc(boundary_definition.entities,parallel=True)
+
+#### Define material properties for the 3 volumes ################
+cubit.cmd('#### DEFINE MATERIAL PROPERTIES #######################')
+
+
+cubit.cmd('block 1 name "elastic 1" ') # elastic material region
+cubit.cmd('block 1 attribute count 6')
+cubit.cmd('block 1 attribute index 1 1') # flag for material: 1 for 1. material
+cubit.cmd('block 1 attribute index 2 2800') # vp
+cubit.cmd('block 1 attribute index 3 1500') # vs
+cubit.cmd('block 1 attribute index 4 2300') # rho
+cubit.cmd('block 1 attribute index 5 150.0') # Qmu
+cubit.cmd('block 1 attribute index 6 0 ') # anisotropy_flag
+
+# optional saves backups
+cubit.cmd('export mesh "top.e" dimension 3 overwrite')
+cubit.cmd('save as "meshing.cub" overwrite')
+# cleanup
+os.system('rm -f topo_brick.stl topo_vol.stl topo_vol2.stl')
+
+#### Export to SESAME format using cubit2specfem3d.py of GEOCUBIT
+
+os.system('mkdir -p MESH')
+cubit2specfem3d.export2SESAME('MESH')
+
+# all files needed by SCOTCH are now in directory MESH
Property changes on: seismo/3D/SPECFEM3D/trunk/examples/Mount_StHelens/mesh_mount.py
___________________________________________________________________
Name: svn:executable
+ *
Added: seismo/3D/SPECFEM3D/trunk/examples/Mount_StHelens/process.sh
===================================================================
--- seismo/3D/SPECFEM3D/trunk/examples/Mount_StHelens/process.sh (rev 0)
+++ seismo/3D/SPECFEM3D/trunk/examples/Mount_StHelens/process.sh 2011-11-30 06:01:19 UTC (rev 19254)
@@ -0,0 +1,82 @@
+#!/bin/bash
+#
+# script runs decomposition,database generation and solver
+# using this example setup
+#
+# prior to running this script, you must create the mesh files
+# in directory MESH/
+#
+
+###################################################
+
+# number of processes
+NPROC=4
+
+##################################################
+
+echo "running example: `date`"
+currentdir=`pwd`
+
+echo
+echo "(will take about 15 minutes)"
+echo
+
+# sets up directory structure in current example directoy
+echo
+echo " setting up example..."
+echo
+
+mkdir -p bin
+mkdir -p in_out_files/OUTPUT_FILES
+mkdir -p in_out_files/DATABASES_MPI
+
+rm -rf in_out_files/OUTPUT_FILES/*
+rm -rf in_out_files/DATABASES_MPI/*
+
+# compiles executables in root directory
+cd ../../
+make >& $currentdir/tmp.log
+cd $currentdir
+
+# links executables
+cd bin/
+rm -f ./x*
+cp ../../../bin/xdecompose_mesh_SCOTCH ./
+cp ../../../bin/xgenerate_databases ./
+cp ../../../bin/xspecfem3D ./
+cd ../
+
+# decomposes mesh
+echo
+echo " decomposing mesh..."
+echo
+./bin/xdecompose_mesh_SCOTCH $NPROC MESH/ in_out_files/DATABASES_MPI/
+
+# stores setup
+cp in_data_files/Par_file in_out_files/OUTPUT_FILES/
+cp in_data_files/CMTSOLUTION in_out_files/OUTPUT_FILES/
+cp in_data_files/STATIONS in_out_files/OUTPUT_FILES/
+
+# runs database generation
+echo
+echo " running database generation..."
+echo
+cd bin/
+mpirun -np $NPROC ./xgenerate_databases
+cd ../
+
+# runs simulation
+echo
+echo " running solver..."
+echo
+cd bin/
+mpirun -np $NPROC ./xspecfem3D
+cd ../
+
+echo
+echo "see results in directory: in_out_files/OUTPUT_FILES/"
+echo
+echo "done"
+echo `date`
+
+
Property changes on: seismo/3D/SPECFEM3D/trunk/examples/Mount_StHelens/process.sh
___________________________________________________________________
Name: svn:executable
+ *
Added: seismo/3D/SPECFEM3D/trunk/examples/Mount_StHelens/ptopo.mean.utm
===================================================================
--- seismo/3D/SPECFEM3D/trunk/examples/Mount_StHelens/ptopo.mean.utm (rev 0)
+++ seismo/3D/SPECFEM3D/trunk/examples/Mount_StHelens/ptopo.mean.utm 2011-11-30 06:01:19 UTC (rev 19254)
@@ -0,0 +1,1156 @@
+553910.515393343 5127401.91810055 1075.98
+554372.604503965 5127406.01729845 986.5
+554834.693588418 5127410.15148573 895.607
+555296.782646479 5127414.32066247 795.018
+555758.871677923 5127418.52482877 838.469
+556220.960682528 5127422.7639847 933.25
+556683.04966007 5127427.03813036 918.393
+557145.138610327 5127431.34726583 793.219
+557684.242350842 5127436.41881391 767.054
+558146.331241213 5127440.80376095 760.232
+558531.405295149 5127444.48461203 785.732
+559070.508937754 5127449.67862542 888.531
+559532.597743476 5127454.16854303 1028.32
+559994.686520535 5127458.69345101 1166.79
+560456.775268707 5127463.25334945 1228.09
+560918.863987769 5127467.84823846 1276.02
+561380.952677498 5127472.47811812 1300.27
+561843.041337669 5127477.14298853 1351.18
+562305.12996806 5127481.84284978 1434.11
+562767.218568446 5127486.57770198 1537.02
+563229.307138605 5127491.34754521 1576.84
+563691.395678312 5127496.15237959 1408.95
+564153.484187343 5127500.9922052 1384.75
+564692.587408812 5127506.68289358 1285
+565077.661112487 5127510.77683053 1168.04
+565539.749528151 5127515.72163046 1248.04
+566078.852639844 5127521.53478925 1360.27
+566540.940986824 5127526.55540453 1277.79
+567003.029301749 5127531.61101166 1445.27
+567465.117584396 5127536.70161077 1417.16
+567927.20583454 5127541.82720194 1256.71
+568389.294051958 5127546.9877853 1105.7
+568851.382236426 5127552.18336093 1328.14
+569313.47038772 5127557.41392896 1509.05
+553916.600170059 5126713.01823802 959.893
+554378.741440846 5126717.11747467 825.796
+554840.882685591 5126721.25170103 756.755
+555303.023904071 5126725.42091719 740
+555765.165096063 5126729.62512324 772.661
+556227.306261345 5126733.86431925 852.286
+556689.447399696 5126738.13850532 777.265
+557228.612026769 5126743.1692794 760.143
+557690.753106005 5126747.51927757 767.306
+558152.894157606 5126751.90426607 848.735
+558538.011679336 5126755.58515196 987.98
+559077.176177011 5126760.77921447 1028.09
+559539.317144372 5126765.26917453 1024.98
+560001.458083207 5126769.7941253 1039.55
+560463.598993295 5126774.35406688 1049.37
+560925.739874413 5126778.94899934 1112.04
+561387.88072634 5126783.57892279 1187.29
+561850.021548852 5126788.24383732 1289.24
+562312.162341727 5126792.94374303 1381.32
+562774.303104742 5126797.67864001 1493.88
+563236.443837676 5126802.44852837 1522.51
+563698.584540305 5126807.25340819 1477.61
+564160.725212407 5126812.09327959 1348.27
+564699.88929098 5126817.78402181 1132.39
+565085.006464139 5126821.87799749 1064.35
+565547.147043324 5126826.8228442 1176.33
+566086.31101265 5126832.63605799 1319.09
+566548.451523482 5126837.65672077 1174.65
+567010.592002415 5126842.71237574 1187.96
+567472.732449224 5126847.80302302 1164.05
+567934.872863688 5126852.92866269 1071.2
+568397.013245583 5126858.08929488 1080.49
+568859.153594688 5126863.28491968 1345.2
+569321.293910778 5126868.51553721 1422.71
+553922.389931585 5126057.45293377 755.786
+554384.580833599 5126061.55220712 686.959
+554846.77170969 5126065.68647049 699.735
+555308.962559637 5126069.85572398 726.408
+555771.15338322 5126074.05996766 743.375
+556233.344180216 5126078.29920163 772.469
+556695.534950404 5126082.57342597 765.327
+557234.757481447 5126087.6042451 772.536
+557696.948192794 5126091.95428222 824.633
+558159.138876633 5126096.33930999 858.388
+558544.297758661 5126100.02022884 844.694
+559083.520160903 5126105.21433786 850.286
+559545.710760891 5126109.70433814 880.898
+560007.901332485 5126114.22932943 905.633
+560470.091875466 5126118.78931184 927.673
+560932.28238961 5126123.38428546 963.357
+561394.472874697 5126128.01425038 1051.31
+561856.663330505 5126132.67920669 1162
+562318.853756814 5126137.3791545 1263.41
+562781.0441534 5126142.1140939 1335.12
+563243.234520044 5126146.88402498 1418.1
+563705.424856523 5126151.68894785 1390.35
+564167.615162616 5126156.52886261 1282.21
+564706.837147691 5126162.21965581 1078.53
+565091.995682757 5126166.31366817 1050.27
+565554.185896362 5126171.25855918 1165.94
+566093.407772694 5126177.07182506 1297.04
+566555.597918262 5126182.09253283 1105.63
+567017.788032076 5126187.1482331 1052.02
+567479.978113917 5126192.23892599 1049.12
+567942.168163561 5126197.3646116 1061.33
+568404.358180786 5126202.52529004 1193.94
+568866.548165372 5126207.7209614 1326.94
+569328.738117097 5126212.95162581 1325.25
+553928.081000191 5125412.99956552 652.804
+554390.320687413 5125417.09887477 674.204
+554852.56034883 5125421.23317436 692.714
+555314.799984224 5125425.40246437 714.592
+555777.039593372 5125429.60674488 723.089
+556239.279176055 5125433.84601599 756.98
+556701.518732054 5125438.12027777 777.878
+557240.798180031 5125443.15114097 775.214
+557703.037677456 5125447.5012162 815.082
+558165.277147499 5125451.8862824 908.327
+558550.476684793 5125455.56723349 1014.67
+559089.756004556 5125460.76138803 1095.71
+559551.99539113 5125465.25142765 1110.22
+560014.234749441 5125469.77645859 1123.92
+560476.474079267 5125474.33648096 1150.51
+560938.71338039 5125478.93149485 1114.52
+561400.952652587 5125483.56150034 1016.39
+561863.191895639 5125488.22649754 1040.27
+562325.431109326 5125492.92648654 1101.5
+562787.670293427 5125497.66146743 1168.53
+563249.909447721 5125502.43144033 1276.45
+563712.148571988 5125507.23640531 1321.18
+564174.387666008 5125512.07636249 1236.75
+564713.666570507 5125517.76720558 1085.86
+565098.865762422 5125521.86125383 1046
+565561.104764376 5125526.80618818 1094.9
+566100.383560628 5125532.61950503 1163.05
+566562.622494854 5125537.64025681 1056.49
+567024.861397473 5125542.69600141 1049
+567487.100268262 5125547.78673893 1049
+567949.339107002 5125552.91246948 1062.08
+568411.577913472 5125558.07319317 1094.12
+568873.816687451 5125563.26891009 1130.71
+569336.055428718 5125568.49962037 1216.18
+553934.262047481 5124712.99061343 706.5
+554396.554720128 5124717.0899615 749.375
+554858.8473671 5124721.22430023 730.929
+555321.139988176 5124725.39362971 729.839
+555783.432583138 5124729.59795004 744.125
+556245.725151767 5124733.83726129 741.554
+556708.017693844 5124738.11156355 765.018
+557170.31020915 5124742.42085691 800.469
+557709.651442878 5124747.49259081 897.875
+558171.943899427 5124751.87769853 969.804
+558557.187592248 5124755.55868449 1081.41
+559096.528729911 5124760.75288822 1184.11
+559558.821103408 5124765.24297036 1134.52
+560021.113448781 5124769.76804417 1166.57
+560483.405765812 5124774.32810973 1221.16
+560945.698054282 5124778.92316714 1207.02
+561407.990313971 5124783.55321649 1153.71
+561870.28254466 5124788.21825788 1191.5
+562332.57474613 5124792.9182914 1228.42
+562794.866918161 5124797.65331716 1252.16
+563257.159060534 5124802.42333524 1276.62
+563719.451173029 5124807.22834575 1371.86
+564181.743255427 5124812.06834878 1219.8
+564721.083979894 5124817.75924578 1060.34
+565106.327329055 5124821.85333283 1050.7
+565568.619319845 5124826.79831404 1068.12
+566107.959936603 5124832.61168597 1063.95
+566570.251860003 5124837.63248533 1049
+567032.543751953 5124842.68827784 1049
+567494.835612232 5124847.77906361 1049
+567957.127440621 5124852.90484274 1068.82
+568419.4192369 5124858.06561533 1208.96
+568881.71100085 5124863.26138151 1196.93
+569344.002732251 5124868.49214136 1254.31
+553940.344341656 5124024.09366793 880.839
+554402.689153195 5124028.193054 803.837
+554865.033939184 5124032.32743106 813
+555327.378699405 5124036.4967992 754.51
+555789.723433642 5124040.7011585 750.696
+556252.068141675 5124044.94050905 763.51
+556714.412823286 5124049.21485094 781.163
+557253.814918148 5124054.24580843 797.589
+557716.159541765 5124058.59596519 821.755
+558178.504138271 5124062.98111357 833.367
+558563.791281158 5124066.66213366 907.02
+559103.193249074 5124071.85638556 987.357
+559565.537762935 5124076.34650934 998.571
+560027.882248812 5124080.87162512 997.122
+560490.226706487 5124085.43173297 994.061
+560952.57113574 5124090.02683299 1071.61
+561414.915536355 5124094.65692529 1102.45
+561877.259908113 5124099.32200994 1097.06
+562339.604250795 5124104.02208706 1163.54
+562801.948564183 5124108.75715673 1193.78
+563264.292848059 5124113.52721906 1156.71
+563726.637102205 5124118.33227414 1229.78
+564188.981326402 5124123.17232208 1193.38
+564728.382883155 5124128.86327188 1050.82
+565113.669684076 5124132.9573969 1050.37
+565576.013817116 5124137.90242399 1050.14
+566115.415266691 5124143.71584986 1049.07
+566577.759332673 5124148.73669581 1049
+567040.103367359 5124153.79253523 1049
+567502.44737053 5124158.88336823 1050.34
+567964.791341968 5124164.00919493 1149.98
+568427.135281455 5124169.17001542 1263.37
+568889.479188771 5124174.36582981 1174.47
+569351.823063699 5124179.59663821 1144.3
+553946.033653806 5123379.64236324 1010.12
+554408.4272355 5123383.74178469 903.755
+554870.820791764 5123387.87619743 834.143
+555333.21432238 5123392.04560155 794.143
+555795.60782713 5123396.24999713 767.911
+556258.001305799 5123400.48938428 789.694
+556720.394758169 5123404.76376306 814.306
+557259.853752404 5123409.79476397 829.5
+557722.247147048 5123414.14495828 826.694
+558184.640514706 5123418.53014451 832.796
+558569.968300318 5123422.21119638 850.061
+559109.427168195 5123427.40549311 878.625
+559571.820453592 5123431.89565566 919.837
+560034.213711134 5123436.4208105 936.571
+560496.606940604 5123440.98095772 949.959
+560959.000141785 5123445.57609742 974.214
+561421.393314459 5123450.20622969 1001.65
+561883.78645841 5123454.87135462 1025.37
+562346.17957342 5123459.57147232 1043.48
+562808.572659272 5123464.30658288 1078.82
+563270.965715749 5123469.0766864 1086.45
+563733.358742633 5123473.88178298 1079.18
+564195.751739706 5123478.72187271 1077.38
+564735.210198326 5123484.41287166 1049.9
+565120.537643553 5123488.50703204 1050.02
+565582.930549892 5123493.45210184 1050.55
+566122.38890183 5123499.26557792 1050.62
+566584.781741421 5123504.28646723 1050.08
+567047.174549861 5123509.34235032 1096.2
+567509.567326932 5123514.43322731 1081.88
+567971.960072417 5123519.55909828 1151.12
+568434.352786099 5123524.71996335 1276.39
+568896.745467759 5123529.91582263 1157.33
+569359.13811718 5123535.14667622 1096.14
+553952.212793017 5122679.63565264 1131.83
+554414.659343785 5122683.73511232 1113.3
+554877.105869251 5122687.86956362 999.982
+555339.552369199 5122692.03900663 853.25
+555801.998843412 5122696.24344144 790.578
+556264.445291676 5122700.48286813 829
+556726.891713773 5122704.75728679 833.125
+557189.338109489 5122709.0666975 860.062
+557728.858870858 5122714.13856953 875.625
+558191.305208672 5122718.52379667 889.768
+558576.677136185 5122722.20488288 887.946
+559116.197802886 5122727.39922808 891.438
+559578.644058855 5122731.88943252 914.518
+560041.09028711 5122736.41462958 947.893
+560503.536487434 5122740.97481935 980.804
+560965.982659613 5122745.57000193 987.219
+561428.428803429 5122750.20017741 1019.14
+561890.874918667 5122754.86534588 1038.7
+562353.321005111 5122759.56550744 1041.81
+562815.767062543 5122764.30066219 1077.3
+563278.213090749 5122769.07081023 1079.5
+563740.659089511 5122773.87595166 1071.91
+564203.105058615 5122778.71608656 1060.14
+564742.625318127 5122784.40713863 1055.5
+565127.996906977 5122788.50133723 1052.11
+565590.442785804 5122793.44645319 1050.54
+566129.962939173 5122799.25998354 1079.55
+566592.408751589 5122804.28091972 1160.32
+567054.854533012 5122809.33685002 1272.48
+567517.300283224 5122814.42777453 1304.53
+567979.74600201 5122819.55369336 1284
+568442.191689153 5122824.71460661 1252.29
+568904.637344436 5122829.91051441 1211.07
+569367.082967643 5122835.14141684 1072.45
+553958.293209192 5121990.74091306 1163.8
+554420.791882758 5121994.84041018 1243.86
+554883.290531148 5121998.97489924 1088.24
+555345.789154148 5122003.14438032 864.184
+555808.287751542 5122007.34885353 786.268
+556270.786323116 5122011.58831893 867.714
+556733.284868655 5122015.86277663 942.98
+557272.866471924 5122020.89387043 887.446
+557735.364960317 5122025.24414506 887.776
+558197.863421994 5122029.62941226 928.306
+558583.278786163 5122033.31053209 947.694
+559122.860264342 5122038.50492473 944.071
+559585.358644583 5122042.9951702 948.224
+560047.856997249 5122047.5204086 967.163
+560510.355322125 5122052.08064003 1011.12
+560972.853618995 5122056.67586459 1020.96
+561435.351887645 5122061.30608237 1043.69
+561897.850127859 5122065.97129346 1073.04
+562360.348339423 5122070.67149797 1090.84
+562822.846522121 5122075.40669599 1102.04
+563285.344675738 5122080.17688761 1101.51
+563747.84280006 5122084.98207295 1089.2
+564210.34089487 5122089.82225208 1084.84
+564749.921967896 5122095.51335616 1079.59
+565135.336995096 5122099.60759218 1066.84
+565597.835000081 5122104.55275333 1050
+566137.415967496 5122110.36633681 1117.32
+566599.913906402 5122115.38731889 1239.73
+567062.41181447 5122120.4432954 1159.16
+567524.909691484 5122125.53426644 1136.68
+567987.407537229 5122130.66023213 1204.33
+568449.905351488 5122135.82119256 1117.53
+568912.403134048 5122141.01714785 1069
+569374.900884691 5122146.24809811 975.339
+553964.078820727 5121335.18048392 1265.79
+554426.627089955 5121339.28001648 1190.18
+554889.175334129 5121343.41454129 1001.47
+555351.723553033 5121347.58405843 943.163
+555814.271746455 5121351.78856799 836.982
+556276.819914179 5121356.02807006 792.653
+556739.368055994 5121360.30256472 943.571
+557279.007521743 5121365.33370203 1048.27
+557741.555606685 5121369.68401428 929.51
+558204.103665039 5121374.0693194 951.02
+558589.56035987 5121377.75047108 975.49
+559129.199701128 5121382.94490865 975.679
+559591.747678435 5121387.43519295 1042.12
+560054.2956283 5121391.96047049 1030.71
+560516.843550507 5121396.52074137 1042.41
+560979.391444842 5121401.11600567 1062.14
+561441.939311092 5121405.74626351 1067.53
+561904.487149043 5121410.41151496 1089.29
+562367.034958481 5121415.11176013 1109.7
+562829.58273919 5121419.84699912 1129.1
+563292.130490958 5121424.61723201 1131.69
+563754.678213571 5121429.42245892 1126.37
+564217.225906813 5121434.26267994 1126.64
+564756.864844856 5121439.95383326 1110.69
+565142.321204331 5121444.0481047 1133.29
+565604.868808177 5121448.99330865 1150.51
+566144.507641114 5121454.80694244 1139.36
+566607.055179199 5121459.82796797 1245.98
+567069.602686592 5121464.88398823 1124.41
+567532.150163079 5121469.97500334 1009.61
+567994.697608447 5121475.10101339 989.898
+568457.24502248 5121480.2620185 990.878
+568919.792404965 5121485.45801876 904.857
+569382.339755686 5121490.6890143 889.589
+553969.765809257 5120690.73190819 1081.86
+554432.362828729 5120694.83147543 1150.29
+554894.959823265 5120698.96603521 1219.92
+555357.556792651 5120703.13558762 1060.63
+555820.153736675 5120707.34013275 912.875
+556282.750655123 5120711.57967068 860.306
+556745.347547784 5120715.8542015 888.408
+557285.043889677 5120720.88538137 1110.75
+557747.640725734 5120725.23573042 978.327
+558210.237535329 5120729.62107265 977.327
+558595.734856291 5120733.30225547 1002.61
+559135.431074282 5120738.49673699 1011.45
+559598.027803214 5120742.98705928 1126.24
+560060.624504832 5120747.51237512 1185.12
+560523.221178923 5120752.07268459 1138.86
+560985.817825275 5120756.66798778 1137.02
+561448.414443675 5120761.2982848 1130.31
+561911.011033909 5120765.96357573 1153.45
+562373.607595764 5120770.66386068 1150.29
+562836.204129027 5120775.39913975 1154.27
+563298.800633486 5120780.16941302 1163.35
+563761.397108926 5120784.9746806 1168.45
+564223.993555135 5120789.81494258 1179.23
+564763.689371817 5120795.50614408 1172.8
+565149.186359007 5120799.60045018 1232.22
+565611.782716243 5120804.54569599 1322.98
+566151.478428315 5120810.359379 1313.46
+566614.074720099 5120815.38044703 1218.02
+567076.670981338 5120820.4365101 1078.8
+567539.267211817 5120825.52756831 926.946
+568001.863411322 5120830.65362177 796.041
+568464.459579642 5120835.81467057 774.02
+568927.055716562 5120841.01071484 742.633
+569389.651821868 5120846.24175467 764.304
+553975.942424266 5119990.72816183 1068.62
+554438.592391181 5119994.82776655 1136.54
+554901.242333288 5119998.96236412 1198.02
+555363.892250375 5120003.13195465 1168.59
+555826.54214223 5120007.33653821 1159.92
+556289.192008643 5120011.5761149 1018.25
+556751.841849401 5120015.8506848 967.661
+557214.491664292 5120020.160248 1185.06
+557754.24974869 5120025.2322995 1042.2
+558216.89950681 5120029.61768183 1025.79
+558602.440951649 5120033.2988983 1069.57
+559142.198943229 5120038.49342732 1082.48
+559604.848621103 5120042.98379068 1160.27
+560067.498271804 5120047.50914789 1218.75
+560530.147895121 5120052.06949906 1202.23
+560992.797490841 5120056.66484428 1232.23
+561455.447058753 5120061.29518365 1221.25
+561918.096598645 5120065.96051725 1260.12
+562380.746110304 5120070.66084519 1265.78
+562843.395593519 5120075.39616757 1230.2
+563306.045048078 5120080.16648447 1254.02
+563768.694473769 5120084.971796 1254.41
+564231.343870379 5120089.81210227 1256.94
+564771.10146272 5120095.50335582 1251.32
+565156.642575509 5120099.59769938 1235.62
+565619.291883606 5120104.54299043 1231.16
+566159.049371877 5120110.35672663 1091.77
+566621.698614859 5120115.37784061 949.286
+567084.347827452 5120120.43394994 826.768
+567546.997009445 5120125.52505473 804.781
+568009.646160624 5120130.65115509 774.804
+568472.295280778 5120135.81225112 683.107
+568934.944369694 5120141.00834294 676.143
+569397.59342716 5120146.23943064 649.734
+553982.02035603 5119301.83633948 1116.21
+554444.722424452 5119305.93598087 1029.41
+554907.424468193 5119310.07061545 1004.8
+555370.126487043 5119314.24024329 1084.98
+555832.828480789 5119318.44486448 1251.38
+556295.530449223 5119322.68447911 1230.31
+556758.232392132 5119326.95908727 1102.12
+557298.051292986 5119331.99035817 1159.89
+557760.753179871 5119336.34078594 1079.9
+558223.455040566 5119340.72620752 1093.82
+558609.03990432 5119344.40745695 1137.35
+559148.85868254 5119349.60203247 1173.68
+559611.560463397 5119354.09243602 1199.88
+560074.262217221 5119358.61783375 1262.47
+560536.9639438 5119363.17822574 1303.63
+560999.665642924 5119367.7736121 1335.25
+561462.367314381 5119372.40399292 1373.92
+561925.068957961 5119377.0693683 1428.59
+562387.770573452 5119381.76973832 1437.43
+562850.472160645 5119386.5051031 1408.51
+563313.173719328 5119391.27546272 1390.18
+563775.87524929 5119396.08081728 1377.86
+564238.576750319 5119400.92116689 1367.46
+564778.39513134 5119406.61247141 1327.8
+565163.979664739 5119410.70685164 1308.86
+565626.681077707 5119415.65218698 1291.8
+566166.499355188 5119421.46597525 1121.89
+566629.200703374 5119426.4871342 949.265
+567091.902021326 5119431.54328882 948.776
+567554.603308834 5119436.63443921 1058.98
+568017.304565685 5119441.76058549 854.714
+568480.00579167 5119446.92172776 787.469
+568942.706986576 5119452.11786612 926.735
+569405.408150192 5119457.34900069 777.964
+553987.705586781 5118657.38982745 908.732
+554450.456390384 5118661.48950299 995.408
+554913.207169424 5118665.624172 1109
+555375.957923692 5118669.79383457 1124.69
+555838.708652978 5118673.99849079 1215.54
+556301.459357072 5118678.23814074 1236.73
+556764.210035765 5118682.5127845 1230.47
+557304.085795189 5118687.54409731 1230.98
+557766.836418127 5118691.89456133 1168.37
+558229.587015001 5118696.28001944 1161.49
+558615.212492333 5118699.96129954 1204.65
+559155.088129713 5118705.15591833 1258.43
+559617.838647132 5118709.6463593 1313.94
+560080.589137647 5118714.17179473 1328.55
+560543.339601049 5118718.73222473 1407.92
+561006.090037127 5118723.32764938 1470.75
+561468.840445671 5118727.95806879 1556.59
+561931.590826472 5118732.62348304 1608.02
+562394.34117932 5118737.32389224 1572.55
+562857.091504005 5118742.05929648 1559.41
+563319.841800316 5118746.82969585 1568.04
+563782.592068044 5118751.63509047 1555.96
+564245.34230698 5118756.47548042 1518.02
+564785.217549067 5118762.16683238 1458.98
+565170.842697631 5118766.26124674 1390.41
+565633.592848927 5118771.2066233 1396.88
+566173.467987972 5118777.02046003 1272.77
+566636.218074797 5118782.04166084 1220.59
+567098.968131534 5118787.09785761 1242.82
+567561.718157972 5118792.18905045 1115.7
+568024.468153901 5118797.31523947 958.184
+568487.218119111 5118802.47642477 1058.45
+568949.968053392 5118807.67260646 1085.29
+569412.717956534 5118812.90378464 881.339
+553993.488267029 5118001.83283976 983.839
+554456.288641175 5118005.93254986 1065.2
+554919.088990878 5118010.06725373 1067.02
+555381.889315931 5118014.23695145 1014.92
+555844.689616125 5118018.44164312 958.357
+556307.489891251 5118022.68132881 962.265
+556770.2901411 5118026.95600862 1004.71
+557310.223733699 5118031.98736385 1050.96
+557773.023928067 5118036.33786455 1141.8
+558235.824096499 5118040.72335964 1173.51
+558621.490883562 5118044.40467079 1181.37
+559161.424354715 5118049.59933339 1255.77
+559624.224444083 5118054.08981222 1364.55
+560087.024506679 5118058.61528582 1448.82
+560549.824542294 5118063.17575428 1536.2
+561012.624550719 5118067.7712177 1657.59
+561475.424531747 5118072.40167615 1794.47
+561938.224485167 5118077.06712976 1751.16
+562401.024410771 5118081.7675786 1662.61
+562863.82430835 5118086.50302279 1644.37
+563326.624177695 5118091.27346241 1770.63
+563789.424018598 5118096.07889756 1799.67
+564252.223830849 5118100.91932835 1705.34
+564792.15690865 5118106.61072833 1621.67
+565177.823368559 5118110.70517722 1525.45
+565640.623093601 5118115.65059552 1424.49
+566180.556068866 5118121.4644813 1359.18
+566643.355729753 5118126.48572448 1381.69
+567106.1553607 5118131.54196392 1279.57
+567568.954961496 5118136.63319972 1158.8
+568031.754531933 5118141.759432 1176.55
+568494.554071802 5118146.92066086 1193.61
+568957.353580894 5118152.1168864 993.245
+569420.153058999 5118157.34810873 854.464
+553999.564360466 5117312.94317525 817.875
+554462.416820366 5117317.04292148 832.125
+554925.26925595 5117321.17766179 793.411
+555388.121667011 5117325.34739626 795.107
+555850.974053342 5117329.55212499 825.906
+556313.826414736 5117333.79184805 872.607
+556776.678750984 5117338.06656554 929.375
+557239.531061879 5117342.37627754 966.484
+557779.52539227 5117347.44850417 1033.75
+558242.377647522 5117351.83403792 1125.7
+558628.087840373 5117355.51538152 1214.84
+559168.082079789 5117360.71008992 1311.23
+559630.934256388 5117365.20060834 1440.21
+560093.786406353 5117369.72612184 1524.23
+560556.638529478 5117374.28663051 1642.75
+561019.490625555 5117378.88213444 1866.3
+561482.342694375 5117383.51263373 2099.5
+561945.194735731 5117388.17812848 2014.2
+562408.046749416 5117392.87861877 1786.41
+562870.898735221 5117397.61410472 1752.43
+563333.750692938 5117402.3845864 1912.07
+563796.60262236 5117407.19006394 2071.27
+564259.454523279 5117412.03053742 1930.39
+564799.448371383 5117417.72198759 1727.45
+565185.158238776 5117421.8164726 1562.71
+565648.010052938 5117426.76193452 1427.96
+566188.003799036 5117432.57587159 1369.02
+566650.855549376 5117437.59715906 1356.89
+567113.707269931 5117442.6534431 1325
+567576.558960492 5117447.74472382 1157.52
+568039.410620851 5117452.87100132 1105.43
+568502.2622508 5117458.03227572 1039.57
+568965.113850131 5117463.22854711 906.179
+569427.965418636 5117468.4598156 753.953
+554005.737803044 5116612.94315271 699.411
+554468.6431832 5116617.04293546 740.51
+554931.548539169 5116621.1777126 839.735
+555394.453870745 5116625.34748422 891.286
+555857.359177723 5116629.5522504 982.571
+556320.264459894 5116633.79201123 1056.22
+556783.169717054 5116638.06676681 1068.04
+557323.225818519 5116643.09821121 1114.77
+557786.131020778 5116647.44878902 1163.02
+558249.036197371 5116651.83436185 1243.63
+558634.790491444 5116655.51573824 1351.41
+559174.846472733 5116660.71049293 1429.27
+559637.75157109 5116665.20105137 1541.02
+560100.656642954 5116669.72660519 1707.18
+560563.56168812 5116674.2871545 1877.53
+561026.466706381 5116678.88269939 2050.41
+561489.37169753 5116683.51323994 2282.35
+561952.27666136 5116688.17877626 2114.29
+562415.181597666 5116692.87930845 1971.05
+562878.086506239 5116697.6148366 1952.98
+563340.991386874 5116702.38536081 2133.16
+563803.896239363 5116707.19088118 2313.08
+564266.8010635 5116712.0313978 2072.98
+564806.856655551 5116717.72289871 1820.37
+565192.610625891 5116721.81742022 1621.41
+565655.51536373 5116726.76292623 1447.63
+566195.570854316 5116732.57691513 1355.88
+566658.475528672 5116737.59824737 1338.35
+567121.380173399 5116742.65457649 1253.86
+567584.284788292 5116747.7459026 1057.18
+568047.189373144 5116752.87222582 884.122
+568510.093927746 5116758.03354623 773.122
+568972.998451893 5116763.22986396 694.551
+569435.902945376 5116768.4611791 658.571
+554011.420708427 5115968.49936982 749.571
+554474.374803836 5115972.59918601 798.122
+554937.328875177 5115976.73399689 860.122
+555400.282922245 5115980.90380252 935.184
+555863.236944834 5115985.108603 1023.16
+556326.19094274 5115989.34839843 1221.61
+556789.144915757 5115993.62318887 1304.67
+557329.257852544 5115998.65467433 1198.32
+557792.211770929 5116003.00528764 1230.49
+558255.165663774 5116007.39089625 1276.67
+558640.960554821 5116011.07230269 1345.88
+559181.073372022 5116016.26709976 1437.93
+559644.027187014 5116020.75769484 1568.55
+560106.980975645 5116025.2832856 1727.53
+560569.934737707 5116029.84387213 1931.9
+561032.888472997 5116034.43945452 2132.55
+561495.842181308 5116039.07003287 2346.69
+561958.795862434 5116043.73560727 2185.22
+562421.74951617 5116048.43617783 2104.79
+562884.70314231 5116053.17174463 2083.35
+563347.656740648 5116057.94230778 2250.41
+563810.610310979 5116062.74786737 2313.9
+564273.563853097 5116067.58842351 2069.96
+564813.676282969 5116073.27997089 1825.04
+565199.47085187 5116077.37452583 1654.04
+565662.424308114 5116082.32007221 1476.8
+566202.536637018 5116088.13410858 1364.84
+566665.490030089 5116093.15548182 1299.27
+567128.443393678 5116098.21185223 1164.37
+567591.396727577 5116103.30321992 1008.64
+568054.350031582 5116108.42958499 859.612
+568517.303305487 5116113.59094756 801.551
+568980.256549085 5116118.78730772 810.408
+569443.209762171 5116124.01866559 979.714
+554017.592887819 5115268.50082952 802.734
+554480.599892659 5115272.60068185 895.232
+554943.606873561 5115276.73552917 909.482
+555406.613830319 5115280.90537156 956.018
+555869.62076273 5115285.11020911 1030
+556332.627670589 5115289.35004191 1139.93
+556795.634553692 5115293.62487004 1326.75
+557258.641411835 5115297.93469359 1202.28
+557798.816047849 5115303.00705152 1227.68
+558261.82285121 5115307.39269879 1264.57
+558647.66183446 5115311.07413769 1339.36
+559187.836380906 5115316.26898057 1416.12
+559650.843106832 5115320.75961524 1520.12
+560113.849806537 5115325.28524591 1664.68
+560576.856479816 5115329.84587265 1815.14
+561039.863126466 5115334.44149557 2006.52
+561502.869746282 5115339.07211475 2253.3
+561965.876339058 5115343.7377303 2394.93
+562428.882904591 5115348.43834231 2365.08
+562891.889442676 5115353.17395088 2327.93
+563354.895953108 5115357.9445561 2334.46
+563817.902435682 5115362.75015809 2110.04
+564280.908890195 5115367.59075692 1890.61
+564821.083051386 5115373.2823545 1669.14
+565206.921714215 5115377.37694557 1528
+565669.928083313 5115382.32253557 1430.96
+566210.102144077 5115388.13662324 1319.16
+566673.108450341 5115393.15804077 1246.23
+567136.11472728 5115398.2144558 1198.12
+567599.120974689 5115403.30586841 1183.53
+568062.127192363 5115408.43227872 1101.64
+568525.133380098 5115413.59368682 1089.48
+568988.139537689 5115418.79009284 1061.29
+569451.14566493 5115424.02149687 1040.48
+554023.666453891 5114579.61413063 817.821
+554486.725522826 5114583.71401833 941.755
+554949.78456795 5114587.84890132 1001
+555412.843589058 5114592.01877969 1057.37
+555875.902585949 5114596.22365352 1210.54
+556338.961558418 5114600.46352289 1261.37
+556802.020506262 5114604.73838791 1233.37
+557342.255914019 5114609.76996114 1168.89
+557805.314807812 5114614.12065035 1212.12
+558268.373676338 5114618.50633547 1211.51
+558654.256047327 5114622.18780613 1277.2
+559194.491336769 5114627.38269384 1367.64
+559657.550128269 5114631.87336728 1447.63
+560120.608893687 5114636.399037 1543.84
+560583.66763282 5114640.9597031 1665.96
+561046.726345464 5114645.55536568 1835.55
+561509.785031416 5114650.18602484 2037.88
+561972.843690472 5114654.85168066 2143.43
+562435.90232243 5114659.55233325 2225.39
+562898.960927084 5114664.28798269 2167.49
+563362.019504232 5114669.0586291 2119.06
+563825.078053671 5114673.86427257 1996.41
+564288.136575196 5114678.7049132 1801.82
+564828.371481424 5114684.39655992 1619.1
+565214.253533691 5114688.49118633 1488.57
+565677.311970254 5114693.43681904 1381.73
+566217.546776588 5114699.25095691 1294.8
+566680.60515065 5114704.27241781 1227.65
+567143.663495542 5114709.3288765 1172.02
+567606.721811061 5114714.42033308 1121
+568069.780097002 5114719.54678765 1116.73
+568532.838353163 5114724.70824034 1173.51
+568995.896579339 5114729.90469124 1084.39
+569458.954775327 5114735.13614046 904.5
+554029.347600173 5113935.17241157 881.446
+554492.455369287 5113939.27233219 919.98
+554955.563114708 5113943.40724837 965.408
+555418.670836234 5113947.57716021 991.306
+555881.778533662 5113951.7820678 1084.29
+556344.886206791 5113956.02197122 1189.16
+556807.993855417 5113960.29687056 1123.18
+557348.28608091 5113965.32848419 1096.91
+557811.393675756 5113969.67920833 1120.71
+558274.501245459 5113974.06492867 1162.98
+558660.42420086 5113977.7464289 1230.73
+559200.716308629 5113982.94135833 1320.54
+559663.823801691 5113987.43206782 1385.33
+560126.931268801 5113991.95777389 1449.43
+560590.038709757 5113996.51847663 1538.96
+561053.146124356 5114001.11417612 1652.36
+561516.253512397 5114005.74487247 1793.78
+561979.360873675 5114010.41056576 1892.37
+562442.46820799 5114015.11125611 1974.21
+562905.575515138 5114019.8469436 1920.82
+563368.682794917 5114024.61762833 1864.92
+563831.790047125 5114029.42331041 1810.82
+564294.897271558 5114034.26398993 1684.75
+564835.188998023 5114039.95568238 1539.41
+565221.111636291 5114044.05034169 1416.43
+565684.218776185 5114048.99601413 1343.92
+566224.510403254 5114054.81019872 1273.34
+566687.617480958 5114059.83169997 1224.27
+567150.724529639 5114064.88819929 1147.71
+567613.831549092 5114069.97969678 1101.54
+568076.938539116 5114075.10619256 1060.78
+568540.045499507 5114080.26768673 1031.29
+569003.152430063 5114085.46417939 1036.06
+569466.259330581 5114090.69567066 976.696
+554035.126125101 5113279.62029953 1063.88
+554498.283429149 5113283.72025345 938.388
+554961.440709624 5113287.85520322 907.51
+555424.597966326 5113292.02514895 911.878
+555887.755199053 5113296.2300907 945.661
+556350.912407605 5113300.47002857 1015.78
+556814.06959178 5113304.74496264 1050.04
+557354.419608905 5113309.77661716 1038.8
+557817.576739573 5113314.12737666 1069.41
+558280.733845227 5113318.51313264 1118.29
+558666.698080686 5113322.19466278 1182.41
+559207.047980687 5113327.38963443 1386.93
+559670.20501009 5113331.88038042 1371.76
+560133.362013674 5113336.40612327 1384.43
+560596.518991237 5113340.96686308 1420.22
+561059.675942578 5113345.56259992 1476.38
+561522.832867495 5113350.19333391 1580.06
+561985.989765788 5113354.85906514 1669.04
+562449.146637253 5113359.55979369 1747.23
+562912.303481691 5113364.29551968 1656.29
+563375.460298898 5113369.0662432 1598.06
+563838.617088675 5113373.87196434 1565.08
+564301.773850819 5113378.71268322 1495.14
+564842.12337146 5113384.40442194 1450.51
+565228.087291403 5113388.49911455 1372.35
+565691.243969439 5113393.44482721 1271.29
+566231.59339119 5113399.25905908 1186.98
+566694.750007353 5113404.28060116 1140.1
+567157.90659464 5113409.3371416 1093.78
+567621.063152849 5113414.4286805 1063.54
+568084.219681778 5113419.55521797 1035.43
+568547.376181226 5113424.71675412 1004.43
+569010.532650991 5113429.91328905 974.837
+569473.68909087 5113435.14482287 967.375
+554041.197851464 5112590.73575862 1383.97
+554504.407203841 5112594.83574735 1186.52
+554967.616532773 5112598.97073225 947.143
+555430.825838061 5112603.14071338 897.429
+555894.035119502 5112607.34569084 919.578
+556357.244376899 5112611.58566472 941.339
+556820.453610049 5112615.8606351 963.036
+557283.662818755 5112620.17060207 989.844
+557824.073531079 5112625.2431288 1014.8
+558287.282686132 5112629.62892203 1041.77
+558673.290296194 5112633.31048345 1114.3
+559213.700920798 5112638.50549922 1222.11
+559676.910000012 5112642.99628337 1211.57
+560140.119053544 5112647.52206466 1252.29
+560603.328081197 5112652.08284321 1288.45
+561066.537082768 5112656.6786191 1338.02
+561529.746058057 5112661.30939243 1390.04
+561992.955006865 5112665.9751633 1457.39
+562456.163928991 5112670.6759318 1475.84
+562919.372824235 5112675.41169802 1402.27
+563382.581692395 5112680.18246208 1331.98
+563845.790533272 5112684.98822406 1317.34
+564308.999346665 5112689.82898407 1324.31
+564849.409593953 5112695.52077117 1301.8
+565235.416890198 5112699.61549857 1257.12
+565698.625619936 5112704.56125327 1242.77
+566239.035768866 5112710.37553455 1159.67
+566702.244437065 5112715.39711931 1086.8
+567165.453076542 5112720.45370273 1047.62
+567628.661687098 5112725.54528492 1016.06
+568091.870268532 5112730.67186597 996.839
+568555.078820643 5112735.83344599 977.661
+569018.287343231 5112741.0300251 953.179
+569481.495836094 5112746.2616034 927.688
+554047.36685604 5111890.74094238 1282.38
+554510.629090642 5111894.8409663 1054.65
+554973.891301928 5111898.97598667 926.061
+555437.153489699 5111903.14600358 897.714
+555900.415653756 5111907.35101713 907.268
+556363.6777939 5111911.59102739 913.408
+556826.939909932 5111915.86603446 928.388
+557367.412347894 5111920.8977749 953.893
+557830.674411 5111925.24860868 969.49
+558293.936449364 5111929.63443956 984.755
+558679.988128957 5111933.31603257 1006.49
+559220.460451067 5111938.51109294 1065.91
+559683.722414009 5111943.00191563 1109.76
+560146.984351412 5111947.52773578 1140.94
+560610.246263076 5111952.08855347 1173.76
+561073.508148803 5111956.68436882 1206.16
+561536.770008393 5111961.3151819 1238.67
+562000.031841647 5111965.98099283 1310.37
+562463.293648366 5111970.68180168 1296.48
+562926.555428351 5111975.41760857 1237.12
+563389.817181401 5111980.18841359 1175.39
+563853.078907318 5111984.99421684 1171.76
+564316.340605903 5111989.83501842 1172.89
+564856.812552774 5111995.52685439 1172.78
+565242.863920276 5111999.62161696 1145.9
+565706.125535666 5112004.56741413 1096.41
+566246.59738472 5112010.38174534 1070.11
+566709.858938909 5112015.40337324 1048.88
+567173.120464535 5112020.46000009 1026.67
+567636.381961398 5112025.551626 994.25
+568099.6434293 5112030.67825109 961.224
+568562.90486804 5112035.83987546 944.306
+569026.166277418 5112041.0364992 926.755
+569489.427657236 5112046.26812244 907.089
+554053.045675208 5111246.30195262 894.821
+554516.356590046 5111250.40200875 1028.43
+554979.667481687 5111254.53706161 1084.63
+555442.978349932 5111258.70711129 947.51
+555906.289194585 5111262.91215787 898.929
+556369.600015447 5111267.15220145 883.898
+556832.910812319 5111271.42724212 903.531
+557373.440044753 5111276.45902209 931.982
+557836.750788969 5111280.80989007 923.102
+558300.06150857 5111285.19575541 948.347
+558686.153755957 5111288.87737736 966.551
+559226.682873129 5111294.07247856 1009.61
+559689.993517691 5111298.56333654 1049.78
+560153.304136845 5111303.08919226 1079.14
+560616.614730392 5111307.65004581 1101.37
+561079.925298134 5111312.24589728 1114.41
+561543.235839872 5111316.87674677 1137.12
+562006.546355408 5111321.54259437 1170.2
+562469.856844544 5111326.24344018 1227.45
+562933.167307082 5111330.9792843 1197.35
+563396.477742823 5111335.75012682 1123.39
+563859.78815157 5111340.55596786 1074.78
+564323.098533122 5111345.39680749 1054.57
+564863.627276967 5111351.08868822 1054.1
+565249.719213854 5111355.18348299 1067.47
+565713.029512636 5111360.12931905 1052.84
+566253.558159162 5111365.94369599 1015.88
+566716.868397055 5111370.96536337 998.347
+567180.17860653 5111376.02202999 984.082
+567643.488787389 5111381.11369595 963.482
+568106.798939434 5111386.24036136 939.918
+568570.109062466 5111391.40202633 923.612
+569033.419156287 5111396.59869095 909.265
+569496.729220697 5111401.83035534 897.804
+554058.821832653 5110590.75261696 755.089
+554522.182262136 5110594.85270568 792.469
+554985.542668544 5110598.98779142 877.755
+555448.903051678 5110603.15787425 843.857
+555912.263411342 5110607.36295427 836.768
+556375.623747339 5110611.60303156 852.898
+556838.984059472 5110615.87810621 861.673
+557379.571059868 5110620.9099262 880.643
+557842.93131962 5110625.26082878 924.204
+558306.291554884 5110629.64672899 950.531
+558692.425065422 5110633.32838022 968.469
+559233.011951159 5110638.52352273 982.607
+559696.372111776 5110643.01441643 994.02
+560159.732247117 5110647.54030815 1010.22
+560623.092356984 5110652.10119797 1019.82
+561086.45244118 5110656.69708599 1023.66
+561549.812499509 5110661.32797232 1057.37
+562013.172531772 5110665.99385703 1071.53
+562476.532537773 5110670.69474024 1120.3
+562939.892517314 5110675.43062203 1171.04
+563403.252470198 5110680.20150251 1139.69
+563866.612396228 5110685.00738178 1011.12
+564329.972295206 5110689.84825993 988.036
+564870.55880956 5110695.54018595 958.041
+565256.692011218 5110699.63501329 1033.22
+565720.051827856 5110704.58088871 1010.92
+566260.638245399 5110710.39531191 968.464
+566723.998001465 5110715.41701926 947.082
+567187.357729262 5110720.47372612 946.082
+567650.717428592 5110725.5654326 922.161
+568114.077099258 5110730.69213881 900.367
+568577.436741062 5110735.85384486 884.898
+569040.796353807 5110741.05055084 879.163
+569504.155937294 5110746.28225688 863.839
+554064.891070887 5109901.87099363 816.5
+554528.303527378 5109905.97111642 819.554
+554991.71596092 5109910.10623651 820.768
+555455.128371318 5109914.27635399 823.071
+555918.540758374 5109918.48146895 824.969
+556381.953121894 5109922.72158148 832.679
+556845.365461681 5109926.99669166 840.875
+557308.77777754 5109931.30679957 848.516
+557849.42544887 5109936.37949219 890.554
+558312.837712212 5109940.76542886 909.339
+558699.014579586 5109944.44711068 993.518
+559239.662165051 5109949.64229637 993.156
+559703.074354157 5109954.1332274 931.339
+560166.486518126 5109958.65915673 924.375
+560629.898656762 5109963.22008447 937.946
+561093.310769868 5109967.81601069 954.438
+561556.722857249 5109972.44693551 963.911
+562020.134918708 5109977.11285902 964.107
+562483.54695405 5109981.8137813 991.859
+562946.958963077 5109986.54970247 1057.02
+563410.370945594 5109991.32062262 1066.18
+563873.782901404 5109996.12654184 972.107
+564337.194830312 5110000.96746024 954.938
+564877.842046441 5110006.65943359 936.893
+565264.018606634 5110010.75429499 893.25
+565727.430453656 5110015.70021153 929.661
+566268.077573506 5110021.51468309 892.844
+566731.489360289 5110026.5364322 904.679
+567194.901118958 5110031.59318111 911.982
+567658.312849317 5110036.68492995 893.344
+568121.724551169 5110041.8116788 862.286
+568585.136224318 5110046.97342778 853.214
+569048.547868568 5110052.17017699 854.286
+569511.959483722 5110057.40192655 845.922
+554070.959671495 5109212.9901178 1002.79
+554534.522863972 5109205.97929926 1045.67
+554997.9881582 5109210.11445407 1041.41
+555461.453429412 5109214.28460656 966.51
+555924.817445925 5109229.60073036 953.464
+556388.383902015 5109222.72990494 882.694
+556851.849103016 5109227.00505101 939.653
+557392.454585698 5109243.1479268 974.393
+557856.023623297 5109236.38793033 975.571
+558319.488748282 5109240.77390382 1016.37
+558705.709667131 5109244.45561656 1104.73
+559246.311680304 5109260.76181527 1096.39
+559709.883975992 5109254.14181462 940.939
+560173.349002164 5109258.66778196 888.327
+560636.814003145 5109263.228748 890.98
+561100.168378061 5109278.9356798 894.893
+561563.743928755 5109272.45567654 880.633
+562027.208852994 5109277.12163923 888.306
+562490.56063344 5109292.93356612 897.411
+562954.138623365 5109286.55856195 887.367
+563417.603469106 5109291.32952217 965.51
+563881.068288291 5109296.13548177 873.98
+564344.416606678 5109312.08740341 845.464
+564885.241971159 5109306.668462 821.122
+565271.462584556 5109310.76335781 831.102
+565734.927295564 5109315.7093159 835.531
+566275.516120109 5109332.63479616 854.839
+566739.100741371 5109326.54562762 871.592
+567202.565364523 5109331.60241903 875.633
+567665.907472147 5109347.80516847 858.786
+568129.49452618 5109341.82100259 849.224
+568592.959064294 5109346.98279494 817.837
+569056.423573671 5109352.17958783 820.429
+569519.762210403 5109368.52233643 812.696
+554076.734037983 5108557.4428817 939.571
+554540.248015315 5108561.54307042 1050.51
+555003.761969947 5108565.678257 1128.69
+555467.275901684 5108569.84844154 1158.08
+555930.789810333 5108574.05362412 1094.27
+556394.3036957 5108578.29380484 1042.71
+556857.817557591 5108582.56898377 1045.41
+557398.583699866 5108587.6009265 970.268
+557862.097510228 5108591.95193521 986.408
+558325.611296499 5108596.33794242 1040.02
+558711.872766518 5108600.01968346 1058.04
+559252.638795998 5108605.21495272 1070.36
+559716.152508837 5108609.70595599 966.857
+560179.666196811 5108614.23195813 851.592
+560643.179859725 5108618.79295923 833.184
+561106.693497386 5108623.3889594 837.232
+561570.2071096 5108628.01995872 805.531
+562033.720696172 5108632.6859573 853.408
+562497.234256908 5108637.38695522 833.625
+562960.747791615 5108642.12295259 763.163
+563424.261300099 5108646.8939495 778.776
+563887.774782164 5108651.69994606 806.388
+564351.288237618 5108656.54094237 801.482
+564892.053901754 5108662.23300731 811.857
+565278.315067912 5108666.32793461 819.531
+565741.828442364 5108671.27393075 818.449
+566282.594011262 5108677.0884959 884.036
+566746.107326127 5108682.11032583 907.653
+567209.620613181 5108687.16715614 891.02
+567673.133872231 5108692.25898694 887.982
+568136.647103083 5108697.38581832 858.306
+568600.160305542 5108702.54765039 810.184
+569063.673479413 5108707.74448327 778.633
+569527.186624503 5108712.97631706 767.464
+554082.801393992 5107868.56346471 775.766
+554546.367382203 5107872.66368692 861.125
+555009.93334784 5107876.79890729 974.857
+555473.49929071 5107880.9691259 1087.61
+555937.065210622 5107885.17434285 1087.72
+556400.631107382 5107889.41455821 1021.86
+556864.196980798 5107893.68977207 958.946
+557327.762830677 5107897.99998453 911.859
+557868.589625529 5107903.07280019 874.464
+558332.15542375 5107907.45884324 948.482
+558718.460237163 5107911.14061438 950.625
+559259.286947555 5107916.33592609 986.453
+559722.852672754 5107920.82696606 885.357
+560186.418373228 5107925.3530052 754.625
+560649.984048782 5107929.91404358 771.625
+561113.549699224 5107934.51008132 755.25
+561577.115324362 5107939.1411185 748.589
+562040.680924001 5107943.80715521 784.536
+562504.24649795 5107948.50819157 765.922
+562967.812046015 5107953.24422765 724.893
+563431.377568004 5107958.01526358 744.804
+563894.943063722 5107962.82129943 802.714
+564358.508532978 5107967.66233531 878.906
+564899.334880072 5107973.3544468 961.696
+565285.639391327 5107977.44940758 874.107
+565749.204780035 5107982.39544417 854.446
+566290.031032425 5107988.21005687 944.438
+566753.596361879 5107993.23192787 1001.57
+567217.161663678 5107998.28879954 1010.32
+567680.72693763 5108003.38067198 970.859
+568144.292183541 5108008.50754529 961.054
+568607.857401219 5108013.66941959 915.5
+569071.422590469 5108018.86629497 865.411
+569534.987751098 5108024.09817155 825.484
+554088.965957237 5107168.57385508 716
+554552.584789611 5107172.67411113 783.755
+555016.203599542 5107176.80936563 909.449
+555479.822386837 5107180.97961866 1038.98
+555943.441151304 5107185.18487032 1064.62
+556407.059892752 5107189.42512068 983.694
+556870.67861099 5107193.70036984 876.163
+557411.567086008 5107198.73239523 819.089
+557875.185753299 5107203.08347541 770.776
+558338.804396772 5107207.46955468 836.122
+558725.153248001 5107211.15135621 835.653
+559266.041611498 5107216.34671082 883.393
+559729.660182368 5107220.83778788 887.898
+560193.278728654 5107225.36386438 743.184
+560656.897250163 5107229.92494043 701.878
+561120.515746704 5107234.52101612 708.929
+561584.134218085 5107239.15209154 674.531
+562047.752664114 5107243.81816679 676.98
+562511.3710846 5107248.51924197 688.964
+562974.98947935 5107253.25531718 679.061
+563438.607848173 5107258.0263925 694.082
+563902.226190876 5107262.83246805 819.163
+564365.844507267 5107267.67354393 982.982
+564906.732509547 5107273.36570243 976.204
+565293.081060348 5107277.46069704 937.184
+565756.699296653 5107282.40677449 965.184
+566297.58720477 5107288.22143522 997.911
+566761.20538216 5107293.24334772 1066.18
+567224.823532054 5107298.30026117 1116.43
+567688.44165426 5107303.39217568 1079.66
+568152.059748585 5107308.51909135 1000.69
+568615.677814839 5107313.6810083 950.531
+569079.295852827 5107318.87792664 915.02
+569542.913862359 5107324.10984646 871.161
+554094.640687122 5106524.13965874 800.393
+554558.308164689 5106528.23994578 944.918
+555021.975619932 5106532.37523152 1069.47
+555485.643052659 5106536.54551607 1098.71
+555949.310462679 5106540.7507995 1051.41
+556412.977849802 5106544.9910819 961.98
+556876.645213839 5106549.26636337 819.571
+557417.590442445 5106554.29842679 719.143
+557881.257755805 5106558.64953986 631
+558344.925045472 5106563.03565228 677.714
+558731.314435295 5106566.71748163 754.245
+559272.259552972 5106571.91287551 808.321
+559735.926770422 5106576.40398652 833.98
+560199.593963418 5106580.93009723 779.102
+560663.261131769 5106585.49120776 685.02
+561126.928275284 5106590.08731819 647.411
+561590.595393773 5106594.71842863 657.735
+562054.262487045 5106599.38453915 626.694
+562517.929554908 5106604.08564987 638.786
+562981.596597172 5106608.82176089 637.265
+563445.263613646 5106613.59287229 678.041
+563908.930604139 5106618.39898418 790.429
+564372.59756846 5106623.24009666 911.179
+564913.542326835 5106628.93229821 933.755
+565299.931417822 5106633.02732379 963.102
+565763.598302481 5106637.97343865 993.286
+566304.542967193 5106643.78814335 1041.23
+566768.209793249 5106648.81009383 1085.78
+567231.876591955 5106653.86704553 1168.22
+567695.543363119 5106658.95899855 1080.55
+568159.210106552 5106664.08595301 997.653
+568622.87682206 5106669.24790901 918.98
+569086.543509454 5106674.44486666 871.449
+569550.210168542 5106679.67682606 875.304
+554100.314858703 5105879.70611659 929.571
+554564.030976679 5105883.80643444 966.633
+555027.747072449 5105887.94175126 1019.86
+555491.463145823 5105892.11206715 900.367
+555955.179196612 5105896.31738218 904.589
+556418.895224627 5105900.55769645 890.673
+556882.611229677 5105904.83301005 766.714
+557423.61320629 5105909.86511129 685.321
+557887.329160934 5105914.21625707 553.551
+558351.045092013 5105918.60240246 607.592
+558737.475016441 5105922.28425949 665.694
+559278.476882717 5105927.47969243 743.429
+559742.192741963 5105931.97083719 735.673
+560205.908576885 5105936.49698194 703.306
+560669.624387294 5105941.05812676 631.837
+561133.340172999 5105945.65427175 623.089
+561597.055933812 5105950.285417 630.653
+562060.771669541 5105954.95156261 588.02
+562524.487379998 5105959.65270868 582.482
+562988.203064993 5105964.38885531 617.959
+563451.918724334 5105969.16000259 761.735
+563915.634357833 5105973.96615063 861.184
+564379.3499653 5105978.80729951 904.821
+564920.351474189 5105984.49954388 920.714
+565306.781101376 5105988.59460026 934.429
+565770.496629605 5105993.54075231 979.694
+566311.498045331 5105999.35550076 1024.16
+566775.213515269 5106004.37748901 1103.57
+567238.928958004 5106009.43447875 1165.96
+567702.644373344 5106014.52647009 1142.18
+568166.359761099 5106019.65346312 984.735
+568630.075121079 5106024.81545795 884.878
+569093.790453095 5106030.0124547 821.347
+569557.505756955 5106035.24445347 776.429
+554106.477550366 5105179.71869995 1002.28
+554570.246496468 5105183.81905107 984.911
+555034.015420493 5105187.95440146 868.929
+555497.784322253 5105192.12475119 744.536
+555961.553201559 5105196.33010036 714.938
+556425.322058223 5105200.57044906 785.964
+556889.090892057 5105204.84579736 720.446
+557352.859702872 5105209.15614536 626.531
+557893.923286143 5105214.22912055 495.786
+558357.692046436 5105218.61530154 492.857
+558744.165995316 5105222.29718847 555.304
+559285.229495986 5105227.49266358 615.938
+559748.998184865 5105231.98384482 602.696
+560212.766849562 5105236.51002631 583.125
+560676.535489889 5105241.07120817 586.875
+561140.304105656 5105245.66739048 593.594
+561604.072696675 5105250.29857334 579.411
+562067.841262757 5105254.96475685 515
+562531.609803714 5105259.66594109 574.719
+562995.378319356 5105264.40212619 706.661
+563459.146809495 5105269.17331222 758.821
+563922.915273942 5105273.97949929 803.946
+564386.683712508 5105278.8206875 846.781
+564927.746857873 5105284.5129781 869.268
+565314.220511242 5105288.60806774 908.125
+565777.988871033 5105293.55425998 961.196
+566319.051923776 5105299.36905567 998.359
+566782.820225615 5105304.39108472 1017.12
+567246.588500408 5105309.44811555 1033.73
+567710.356747966 5105314.54014825 1028.55
+568174.1249681 5105319.66718294 931.661
+568637.893160621 5105324.82921972 838.018
+569101.661325339 5105330.0262587 771.75
+569565.429462067 5105335.25829998 734.875
Added: seismo/3D/SPECFEM3D/trunk/examples/Mount_StHelens/read_topo.py
===================================================================
--- seismo/3D/SPECFEM3D/trunk/examples/Mount_StHelens/read_topo.py (rev 0)
+++ seismo/3D/SPECFEM3D/trunk/examples/Mount_StHelens/read_topo.py 2011-11-30 06:01:19 UTC (rev 19254)
@@ -0,0 +1,68 @@
+#!/usr/bin/env python
+
+import cubit
+
+import os
+import sys
+import fileinput
+import string
+import math
+
+print sys.path
+
+#############################################################
+# USER PARAMETERS
+
+# topography file, data points per lon-increment
+inputFile = 'ptopo.mean.utm'
+
+# X coordinate in topography file repeats after line
+nstep = 34
+
+#############################################################
+
+# converts xyz to utm
+os.system('./convert_lonlat2utm.pl ptopo.mean.xyz 10 > ptopo.mean.utm ')
+
+# reads in points
+print '#reading from file: ',inputFile
+
+count = 0
+cubit.cmd('reset')
+cubit.cmd('echo off')
+# creates point vertices
+for line in fileinput.input( inputFile ):
+ count = count + 1
+ #print '# '+str(count)+' point: '+line
+ lineitems = line.split()
+ x = lineitems[0]
+ y = lineitems[1]
+ z = lineitems[2]
+ xyz = str(x) + ' ' + str(y) + ' ' + str(z)
+ #print '#point: ',xyz
+ cubit.cmd('create vertex '+ xyz )
+fileinput.close()
+
+# creates smooth spline curves for surface
+for i in range(1,count+1):
+ if i > 1 :
+ if i % nstep == 0 :
+ cubit.cmd('create curve spline vertex '+str(i-nstep+1) + ' to ' + str(i) )
+cubit.cmd('echo on')
+
+# creates surface
+cubit.cmd('create surface skin curve all')
+
+# cleans up
+cubit.cmd('merge all ')
+cubit.cmd('delete vertex all')
+cubit.cmd('delete curve all')
+
+print '#done points: '+str(count)
+
+# saves and exports surface
+# cubit file (uses ACIS file format)
+cubit.cmd('save as "topo.cub" overwrite')
+# export surface as STL file
+cubit.cmd('export stl ascii "topo.stl" overwrite')
+
Property changes on: seismo/3D/SPECFEM3D/trunk/examples/Mount_StHelens/read_topo.py
___________________________________________________________________
Name: svn:executable
+ *
Added: seismo/3D/SPECFEM3D/trunk/examples/Mount_StHelens/run_boundary_definition.py
===================================================================
--- seismo/3D/SPECFEM3D/trunk/examples/Mount_StHelens/run_boundary_definition.py (rev 0)
+++ seismo/3D/SPECFEM3D/trunk/examples/Mount_StHelens/run_boundary_definition.py 2011-11-30 06:01:19 UTC (rev 19254)
@@ -0,0 +1,18 @@
+#!python
+#!/usr/bin/env python
+
+import cubit
+import boundary_definition
+import cubit2specfem3d
+
+import os
+import sys
+
+
+###### This is boundary_definition.py of GEOCUBIT
+#..... which extracts the bounding faces and defines them into blocks
+reload(boundary_definition)
+boundary_definition.entities=['face']
+boundary_definition.define_bc(boundary_definition.entities,parallel=True)
+
+
Property changes on: seismo/3D/SPECFEM3D/trunk/examples/Mount_StHelens/run_boundary_definition.py
___________________________________________________________________
Name: svn:executable
+ *
Added: seismo/3D/SPECFEM3D/trunk/examples/Mount_StHelens/run_cubit2specfem3d.py
===================================================================
--- seismo/3D/SPECFEM3D/trunk/examples/Mount_StHelens/run_cubit2specfem3d.py (rev 0)
+++ seismo/3D/SPECFEM3D/trunk/examples/Mount_StHelens/run_cubit2specfem3d.py 2011-11-30 06:01:19 UTC (rev 19254)
@@ -0,0 +1,26 @@
+#!python
+#!/usr/bin/env python
+
+import cubit
+import boundary_definition
+import cubit2specfem3d
+
+import os
+import sys
+
+
+###### This is boundary_definition.py of GEOCUBIT
+#..... which extracts the bounding faces and defines them into blocks
+#reload(boundary_definition)
+#boundary_definition.entities=['face']
+#boundary_definition.define_bc(boundary_definition.entities,parallel=True)
+
+
+#### 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
+
Property changes on: seismo/3D/SPECFEM3D/trunk/examples/Mount_StHelens/run_cubit2specfem3d.py
___________________________________________________________________
Name: svn:executable
+ *
Modified: seismo/3D/SPECFEM3D/trunk/examples/homogeneous_halfspace/process.sh
===================================================================
--- seismo/3D/SPECFEM3D/trunk/examples/homogeneous_halfspace/process.sh 2011-11-30 00:51:28 UTC (rev 19253)
+++ seismo/3D/SPECFEM3D/trunk/examples/homogeneous_halfspace/process.sh 2011-11-30 06:01:19 UTC (rev 19254)
@@ -40,9 +40,10 @@
# links executables
cd bin/
-ln -s ../../../bin/xdecompose_mesh_SCOTCH
-ln -s ../../../bin/xgenerate_databases
-ln -s ../../../bin/xspecfem3D
+rm -f ./x*
+cp ../../../bin/xdecompose_mesh_SCOTCH ./
+cp ../../../bin/xgenerate_databases
+cp ../../../bin/xspecfem3D
cd ../
# decomposes mesh
Modified: seismo/3D/SPECFEM3D/trunk/examples/layered_halfspace/process.sh
===================================================================
--- seismo/3D/SPECFEM3D/trunk/examples/layered_halfspace/process.sh 2011-11-30 00:51:28 UTC (rev 19253)
+++ seismo/3D/SPECFEM3D/trunk/examples/layered_halfspace/process.sh 2011-11-30 06:01:19 UTC (rev 19254)
@@ -54,9 +54,10 @@
# links executables
cd bin/
-ln -s ../../../bin/xdecompose_mesh_SCOTCH
-ln -s ../../../bin/xgenerate_databases
-ln -s ../../../bin/xspecfem3D
+rm -f ./x*
+cp ../../../bin/xdecompose_mesh_SCOTCH ./
+cp ../../../bin/xgenerate_databases ./
+cp ../../../bin/xspecfem3D ./
cd ../
# decomposes mesh
Modified: seismo/3D/SPECFEM3D/trunk/examples/tomographic_model/process.sh
===================================================================
--- seismo/3D/SPECFEM3D/trunk/examples/tomographic_model/process.sh 2011-11-30 00:51:28 UTC (rev 19253)
+++ seismo/3D/SPECFEM3D/trunk/examples/tomographic_model/process.sh 2011-11-30 06:01:19 UTC (rev 19254)
@@ -45,9 +45,10 @@
# links executables
cd bin/
-ln -s ../../../bin/xdecompose_mesh_SCOTCH
-ln -s ../../../bin/xgenerate_databases
-ln -s ../../../bin/xspecfem3D
+rm -f ./x*
+cp ../../../bin/xdecompose_mesh_SCOTCH ./
+cp ../../../bin/xgenerate_databases ./
+cp ../../../bin/xspecfem3D ./
cd ../
# decomposes mesh
Modified: seismo/3D/SPECFEM3D/trunk/examples/waterlayered_halfspace/process.sh
===================================================================
--- seismo/3D/SPECFEM3D/trunk/examples/waterlayered_halfspace/process.sh 2011-11-30 00:51:28 UTC (rev 19253)
+++ seismo/3D/SPECFEM3D/trunk/examples/waterlayered_halfspace/process.sh 2011-11-30 06:01:19 UTC (rev 19254)
@@ -40,9 +40,10 @@
# links executables
cd bin/
-ln -s ../../../bin/xdecompose_mesh_SCOTCH
-ln -s ../../../bin/xgenerate_databases
-ln -s ../../../bin/xspecfem3D
+rm -f ./x*
+cp ../../../bin/xdecompose_mesh_SCOTCH ./
+cp ../../../bin/xgenerate_databases ./
+cp ../../../bin/xspecfem3D ./
cd ../
# decomposes mesh
Modified: seismo/3D/SPECFEM3D/trunk/src/shared/create_movie_shakemap_AVS_DX_GMT.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/shared/create_movie_shakemap_AVS_DX_GMT.f90 2011-11-30 00:51:28 UTC (rev 19253)
+++ seismo/3D/SPECFEM3D/trunk/src/shared/create_movie_shakemap_AVS_DX_GMT.f90 2011-11-30 06:01:19 UTC (rev 19254)
@@ -104,8 +104,8 @@
!!!! NL NL for external meshes
!--------------------------------------------
! muting source region
+ logical, parameter :: MUTE_SOURCE = .false.
real(kind=CUSTOM_REAL), parameter :: RADIUS_TO_MUTE = 1000._CUSTOM_REAL
- logical, parameter :: MUTE_SOURCE = .true.
real(kind=CUSTOM_REAL), parameter :: X_SOURCE_EXT_MESH = -9023.021484375
real(kind=CUSTOM_REAL), parameter :: Y_SOURCE_EXT_MESH = 6123.611328125
real(kind=CUSTOM_REAL), parameter :: Z_SOURCE_EXT_MESH = 17.96331405639648
@@ -684,7 +684,7 @@
if(.not. mask_point(ibool_number)) then
call utm_geo(long,lat,xp_save(ilocnum+ieoff),yp_save(ilocnum+ieoff), &
UTM_PROJECTION_ZONE,IUTM2LONGLAT,SUPPRESS_UTM_PROJECTION)
- write(11,*) long,lat,field_display(ilocnum+ieoff)
+ write(11,*) sngl(long),sngl(lat),sngl(field_display(ilocnum+ieoff))
endif
mask_point(ibool_number) = .true.
enddo
@@ -704,10 +704,12 @@
ipoin = ipoin + 1
ireorder(ibool_number) = ipoin
if(USE_OPENDX) then
- write(11,*) xp_save(ilocnum+ieoff),yp_save(ilocnum+ieoff),zp_save(ilocnum+ieoff)
+ write(11,*) sngl(xp_save(ilocnum+ieoff)),sngl(yp_save(ilocnum+ieoff)),sngl(zp_save(ilocnum+ieoff))
else if(USE_AVS) then
- write(11,'(i9,3f16.6)') ireorder(ibool_number),xp_save(ilocnum+ieoff), &
- yp_save(ilocnum+ieoff),zp_save(ilocnum+ieoff)
+ !write(11,'(i9,3f16.6)') ireorder(ibool_number),xp_save(ilocnum+ieoff), &
+ ! yp_save(ilocnum+ieoff),zp_save(ilocnum+ieoff)
+ write(11,*) ireorder(ibool_number),sngl(xp_save(ilocnum+ieoff)), &
+ sngl(yp_save(ilocnum+ieoff)),sngl(zp_save(ilocnum+ieoff))
endif
endif
mask_point(ibool_number) = .true.
@@ -757,13 +759,15 @@
if(plot_shaking_map) then
write(11,*) sngl(field_display(ilocnum+ieoff))
else
- write(11,"(f7.2)") field_display(ilocnum+ieoff)
+ !write(11,"(f7.2)") field_display(ilocnum+ieoff)
+ write(11,*) sngl(field_display(ilocnum+ieoff))
endif
else
if(plot_shaking_map) then
- write(11,*) ireorder(ibool_number),field_display(ilocnum+ieoff)
+ write(11,*) ireorder(ibool_number),sngl(field_display(ilocnum+ieoff))
else
- write(11,"(i10,1x,f7.2)") ireorder(ibool_number),field_display(ilocnum+ieoff)
+ !write(11,"(i10,1x,f7.2)") ireorder(ibool_number),field_display(ilocnum+ieoff)
+ write(11,*) ireorder(ibool_number),sngl(field_display(ilocnum+ieoff))
endif
endif
endif
More information about the CIG-COMMITS
mailing list