Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

Fixed orientation issue / rake is required to be set in set_corners() #323

Merged
merged 5 commits into from
Aug 12, 2018
Merged

Fixed orientation issue / rake is required to be set in set_corners() #323

Changes from all commits
Commits
Show all changes
5 commits
Select commit Hold shift + click to select a range
a32a328
reverse orientation when normal vector points downwards
dsrim Aug 2, 2018
be6e92f
- use when orientation of ``self._corners`` is not counter-clockwise,…
Aug 3, 2018
8661e17
Merge branch 'master' of https://github.com/clawpack/geoclaw into tri…
Aug 7, 2018
316b6cb
``self.rake`` is no longer required to be set in ``self.set_corners``…
Aug 7, 2018
bc87936
removed debug statements and comments
Aug 9, 2018
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
45 changes: 24 additions & 21 deletions src/python/geoclaw/dtopotools.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -1316,6 +1316,7 @@ def __init__(self):
self._corners = None
self._gauss_pts = None
self.n_gauss_pts = 4
self._fix_orientation = False


def convert_to_standard_units(self, input_units, verbose=False):
Expand Down Expand Up @@ -1523,17 +1524,17 @@ def calculate_geometry(self):
self._centers[0][1]
- up_strike[1])

def calculate_geometry_triangles(self,rake=90.):
def calculate_geometry_triangles(self):
r"""
Calculate geometry for triangular subfaults

- Uses *corners* to calculate *centers*, *longitude*, *latitude*,
*depth*, *strike*, *dip*, *rake*, *length*, *width*.
*depth*, *strike*, *dip*, *length*, *width*.

- sets *coordinate_specification* as "triangular"

- Note that calculate_geometry() computes
long/lat/strik/dip/rake/length/width to calculate centers/corners
long/lat/strike/dip/length/width to calculate centers/corners

"""

Expand All @@ -1552,9 +1553,8 @@ def calculate_geometry_triangles(self,rake=90.):

x0 = numpy.array(self.corners)
x = x0.copy()
x = numpy.array(x0)

x[:,2] = - numpy.abs(x[:,2]) # set depth to be always negative(lazy)
x[:,2] = - numpy.abs(x[:,2]) # set depth to be always negative

if 0:
# old coordinate transform
Expand All @@ -1581,6 +1581,7 @@ def calculate_geometry_triangles(self,rake=90.):
normal = cross(v1,v2)
if normal[2] < 0:
normal = -normal
self._fix_orientation = True
strikev = cross(normal,e3) # vector in strike direction

a = normal[0]
Expand All @@ -1589,19 +1590,16 @@ def calculate_geometry_triangles(self,rake=90.):

#Compute strike
strike_deg = rad2deg(numpy.arctan(-b/a))
print('+++ initial strike_deg = %g' % strike_deg)

#Compute dip
beta = deg2rad(strike_deg + 90)
m = numpy.array([sin(beta),cos(beta),0]) #Points in dip direction
n = numpy.array([a,b,c]) #Normal to the plane
#dip_deg = abs(rad2deg(numpy.arcsin(m.dot(n)/(norm(m)*norm(n)))))

if abs(c) < 1e-8:
dip_deg = 90. # vertical fault
else:
dip_deg = rad2deg(numpy.arcsin(m.dot(n)/(norm(m)*norm(n))))
print('+++ initial dip_deg = %g' % dip_deg)

# dip should be between 0 and 90. If negative, reverse strike:
if dip_deg < 0:
Expand All @@ -1616,7 +1614,6 @@ def calculate_geometry_triangles(self,rake=90.):

self.strike = strike_deg
self.dip = dip_deg
self.rake = rake # set default rake to 90 degrees

# find the center line
xx = numpy.zeros((3,3))
Expand Down Expand Up @@ -1649,18 +1646,22 @@ def calculate_geometry_triangles(self,rake=90.):
raise ValueError("Invalid coordinate specification %s." \
% self.coordinate_specification)

def set_corners(self,corners,rake=90.,projection_zone=None):
def set_corners(self,corners,projection_zone=None):
r"""
set three corners for a triangular fault.
Input *corners* should be iterable of length 3.
Set three corners for a triangular fault.
:Inputs:

- *corners* should be iterable of length 3.
- *rake* should be between -180. and 180.

"""

if len(corners) == 3:
self._corners =\
[corners[0],corners[1],corners[2]]
self._projection_zone = projection_zone
self.coordinate_specification = 'triangular'
self.calculate_geometry_triangles(rake=rake)
self.calculate_geometry_triangles()
else:
raise ValueError("Expected input of length 3")

Expand Down Expand Up @@ -1873,9 +1874,13 @@ def okada(self, x, y):
Odepth = abs(O2_list[k][2])

if reverse_list[k]:
sgn = (-1.)**(floor(j/3)+1)
sgn = (-1.)**(floor(j/3))
else:
sgn = (-1.)**floor(j/3)
sgn = (-1.)**floor(j/3+1)

# fix orientation
if self._fix_orientation:
sgn *= -1.

Y1,Y2,Y3,Z1,Z2,Z3,Yb1,Yb2,Yb3,Zb1,Zb2,Zb3 = \
self._get_halfspace_coords(X1,X2,X3,alpha,beta,Olong,Olat,Odepth)
Expand Down Expand Up @@ -1906,7 +1911,7 @@ def okada(self, x, y):
dZ = -v31*burgersv[0] - v32*burgersv[1] + v33*burgersv[2]

dtopo = DTopography()
dtopo.X = X1 # DR: X1, X2 varname confusing?
dtopo.X = X1
dtopo.Y = X2
dtopo.dX = numpy.array(dX, ndmin=3)
dtopo.dY = numpy.array(dY, ndmin=3)
Expand Down Expand Up @@ -1940,7 +1945,7 @@ def _get_leg_angles(self):
y[:,1] = LAT2METER * x[:,1]
y[:,2] = - numpy.abs(x[:,2]) # force sign

v_list = [y[1,:] - y[0,:], y[2,:] - y[1,:], y[0,:] - y[2,:]]
v_list = [y[0,:] - y[1,:], y[1,:] - y[2,:], y[2,:] - y[0,:]]

e3 = numpy.array([0.,0.,-1.])

Expand All @@ -1961,10 +1966,8 @@ def _get_leg_angles(self):
l = j
reverse_list[j] = True

O1 = x[k,:].tolist() # set origin for the vector v
O2 = x[l,:].tolist() # set dest. for the vector v
O1_list.append(O1)
O2_list.append(O2)
O1_list.append(x[k,:].copy()) # set origin for the vector v
O2_list.append(x[l,:].copy()) # set dest. for the vector v

alpha = numpy.arctan2(vn[0],vn[1])
alpha_list.append(alpha)
Expand Down