From 8197f425188a17a2cf0d7639e344eabc48e789d0 Mon Sep 17 00:00:00 2001
From: Surendra Somala <surendra@caltech.edu>
Date: Thu, 23 May 2013 19:19:41 +0000
Subject: [PATCH] Added post- and pre- processing utilities for faults with
 split nodes

---
 utils/FSEM3D_snapshot.m   | 61 +++++++++++++++++++++++++++++++++++++++
 utils/critical_timestep.m | 31 ++++++++++++++++++++
 2 files changed, 92 insertions(+)
 create mode 100644 utils/FSEM3D_snapshot.m
 create mode 100644 utils/critical_timestep.m

diff --git a/utils/FSEM3D_snapshot.m b/utils/FSEM3D_snapshot.m
new file mode 100644
index 000000000..411ef0a3c
--- /dev/null
+++ b/utils/FSEM3D_snapshot.m
@@ -0,0 +1,61 @@
+%FSEM3D_SNAPSHOT reads fault data at a given time
+%
+% d = FSEM3D_snapshot(isnap, [dir, fault])
+%
+% INPUTS	isnap	snapshot index, as in Snapshot*.bin file names
+%		dir	["."] directory containing the SPECFEM3D output data Snapshot*.bin
+%		fault	[1] fault id
+%
+% OUTPUTS	d	structure containing fault fields:
+%			X,Y,Z	  fault node coordinates (km)
+%			Dx,Dz	  slip (m)
+%			Vx,Vz	  slip rate (m/s)
+%			Tx,Ty,Tz  stress change (MPa)
+%			S	  slip "state" variable in friction law (m)
+%			Sg	  strength relative to initial stress (MPa)
+%			Trup	  rupture time (s)
+%			Tpz	  process zone time, when slip=Dc (s)
+%
+% Jean-Paul Ampuero	ampuero@erdw.ethz.ch modified by 
+% Percy Galvez  percy.galvez@sed.ethz.ch 19/01/2011.
+%
+% WARNING : Works only for single precision snapshot files.
+
+function d = FSEM3D_snapshot(isnap,DATA_DIR,fault)
+
+NDAT = 14; 
+
+if nargin<2, DATA_DIR = '.'; end
+if nargin<3, fault = 1; end
+
+BinFile = sprintf('%s/Snapshot%u_F%u.bin',DATA_DIR,isnap,fault);
+
+if ~exist(BinFile,'file'), error(sprintf('File %s does not exist',BinFile)), end
+fid=fopen(BinFile);
+BinRead = fread(fid,[1,inf],'single')' ;
+fclose(fid);
+
+BinRead = reshape( BinRead(:),length(BinRead)/NDAT, NDAT);
+BinRead = BinRead(2:end-1,:);
+
+% Reorder fault nodes (lexicographic z,x)
+%[LOC,IND] = sortrows( BinRead(:,[1 3]) );
+%BinRead = BinRead(IND,:);
+
+d.X  = BinRead(:,1)/1e3; % in km
+d.Y  = BinRead(:,2)/1e3; % in km
+d.Z  = BinRead(:,3)/1e3; % in km
+d.Dx = BinRead(:,4);
+d.Dz = BinRead(:,5);
+d.Vx = BinRead(:,6);
+d.Vz = BinRead(:,7);
+d.Tx = BinRead(:,8); % in MPa
+d.Ty = BinRead(:,9);
+d.Tz = BinRead(:,10); % in MPa
+d.S  = BinRead(:,11);
+d.Sg = BinRead(:,12); % in MPa
+d.Trup = BinRead(:,13);
+d.Tpz = BinRead(:,14);
+
+
+return
diff --git a/utils/critical_timestep.m b/utils/critical_timestep.m
new file mode 100644
index 000000000..cf4b979f0
--- /dev/null
+++ b/utils/critical_timestep.m
@@ -0,0 +1,31 @@
+% CRITICAL_TIMESTEP	Computes the critical timestep in 3D, assuming:
+%                       purely elastic medium
+%                       the leapfrog time scheme
+%                       a cube element (same size for all faces)
+%
+% SYNTAX	dtc = critical_timestep(csp,h,ngll)
+%
+% INPUT		cp	P wave speed (in m/s)
+%		h	element size (in m)
+%		ngll	number of GLL nodes per spectral element edge (2 to 20)
+%
+% OUTPUT	dtc	critical timestep
+%
+function dtc = critical_timestep(csp,h,ngll)
+
+DIM = 3; % dimension
+
+if ngll<2 | ngll>20, error('ngll must be from 2 to 20'); end
+
+% tabulated values of critical frequency (non-dimensional, 1D)
+% Omega_max(p) = Omega_max(ngll-1)
+Omega_max = [2.0000000e+00 4.8989795e+00 8.6203822e+00 1.3540623e+01 1.9797952e+01 2.7378050e+01 ...
+     3.6256848e+01 4.6421894e+01 5.7867306e+01 7.0590158e+01 8.4588883e+01 9.9862585e+01 ...
+     1.1641072e+02 1.3423295e+02 1.5332903e+02 1.7369883e+02 1.9534221e+02 2.1825912e+02 2.4244948e+02];
+
+% stability factor por leapfrog timescheme
+C = 2;
+
+% critical time step, 
+% assumes a cube element
+dtc = C*h/csp/sqrt(DIM)/Omega_max(ngll-1);