| 1 |
######################################################## |
| 2 |
# |
| 3 |
# Copyright (c) 2003-2009 by University of Queensland |
| 4 |
# Earth Systems Science Computational Center (ESSCC) |
| 5 |
# http://www.uq.edu.au/esscc |
| 6 |
# |
| 7 |
# Primary Business: Queensland, Australia |
| 8 |
# Licensed under the Open Software License version 3.0 |
| 9 |
# http://www.opensource.org/licenses/osl-3.0.php |
| 10 |
# |
| 11 |
######################################################## |
| 12 |
|
| 13 |
__copyright__="""Copyright (c) 2003-2009 by University of Queensland |
| 14 |
Earth Systems Science Computational Center (ESSCC) |
| 15 |
http://www.uq.edu.au/esscc |
| 16 |
Primary Business: Queensland, Australia""" |
| 17 |
__license__="""Licensed under the Open Software License version 3.0 |
| 18 |
http://www.opensource.org/licenses/osl-3.0.php""" |
| 19 |
__url__="https://launchpad.net/escript-finley" |
| 20 |
|
| 21 |
from esys.escript import * |
| 22 |
import numpy |
| 23 |
import math |
| 24 |
|
| 25 |
class FaultSystem: |
| 26 |
""" |
| 27 |
The FaultSystem class defines a system of faults in the Earth crust. |
| 28 |
|
| 29 |
A fault system is defined by set of faults index by a tag. Each fault is defined as a polygon describing the top of |
| 30 |
the fault and - in case of a 3D model - a polygon describing the bottom of the fault. This class provides a mechanism |
| 31 |
to parametrise each fault with the domain [0,w0_max] x [0, w1_max] (to [0,w0_max] in the 2D case). |
| 32 |
""" |
| 33 |
NOTAG="__NOTAG__" |
| 34 |
def __init__(self,dim=3): |
| 35 |
""" |
| 36 |
Sets up the fault system |
| 37 |
|
| 38 |
:param dim: spatial dimension |
| 39 |
:type dim: ``int`` of value 2 or 3 |
| 40 |
""" |
| 41 |
if not (dim == 2 or dim == 3): |
| 42 |
raise ValueError,"only dimension2 2 and 3 are supported." |
| 43 |
self.__dim=dim |
| 44 |
self.__faults={} |
| 45 |
self.__length={} |
| 46 |
self.__depth={} |
| 47 |
self.__offsets={} |
| 48 |
self.__w1_max={} |
| 49 |
self.__w0_max={} |
| 50 |
self.__center=None |
| 51 |
self.__orientation = None |
| 52 |
|
| 53 |
def getDim(self): |
| 54 |
""" |
| 55 |
returns the spatial dimension |
| 56 |
:rtype: ``int`` |
| 57 |
""" |
| 58 |
return self.__dim |
| 59 |
def getLength(self,tag=None): |
| 60 |
""" |
| 61 |
returns the unrolled length of fault ``tag`` |
| 62 |
:rtype: ``float`` |
| 63 |
""" |
| 64 |
if tag==None: tag=self.NOTAG |
| 65 |
return self.__length[tag] |
| 66 |
|
| 67 |
def getDepth(self,tag=None): |
| 68 |
""" |
| 69 |
returns the medium depth of fault ``tag`` |
| 70 |
:rtype: ``float`` |
| 71 |
""" |
| 72 |
if tag==None: tag=self.NOTAG |
| 73 |
return self.__depth[tag] |
| 74 |
|
| 75 |
def getTags(self): |
| 76 |
""" |
| 77 |
returns a list of the tags used by the fault system |
| 78 |
:rtype: ``list`` |
| 79 |
""" |
| 80 |
return self.__faults.keys() |
| 81 |
|
| 82 |
def getW0Range(self,tag=None): |
| 83 |
""" |
| 84 |
returns the range of the parameterization in ``w0`` |
| 85 |
:rtype: two ``float`` |
| 86 |
""" |
| 87 |
return self.getW0Offsets(tag)[0], self.getW0Offsets(tag)[-1] |
| 88 |
|
| 89 |
def getW1Range(self,tag=None): |
| 90 |
""" |
| 91 |
returns the range of the parameterization in ``w1`` |
| 92 |
:rtype: two ``float`` |
| 93 |
""" |
| 94 |
if tag==None: tag=self.NOTAG |
| 95 |
return -self.__w1_max[tag],0 |
| 96 |
|
| 97 |
def getW0Offsets(self, tag=None): |
| 98 |
""" |
| 99 |
returns the offsets for the parametrization of fault ``tag``. |
| 100 |
|
| 101 |
:return: the offsets in the parametrization |
| 102 |
:rtype: ``list`` of ``float`` |
| 103 |
""" |
| 104 |
if tag==None: tag=self.NOTAG |
| 105 |
return self.__offsets[tag] |
| 106 |
|
| 107 |
|
| 108 |
def getFaultSegments(self, tag=None): |
| 109 |
""" |
| 110 |
returns the polygons used to describe fault tagged by ``tag`` |
| 111 |
|
| 112 |
:param tag: the tag of the fault |
| 113 |
:type tag: ``float`` or ``str`` |
| 114 |
:return: the list of vertices defining the top and the list of the vertices defining the bottom of the fault indexed by ``tag``. The coordinates are `numpy.ndarray'. |
| 115 |
""" |
| 116 |
if tag==None: tag=self.NOTAG |
| 117 |
return self.__faults[tag] |
| 118 |
|
| 119 |
def getCenterOnSurface(self): |
| 120 |
""" |
| 121 |
returns the center point of the fault system at the surface |
| 122 |
:rtype: `numpy.ndarray` |
| 123 |
""" |
| 124 |
if self.__center == None: |
| 125 |
self.__center=numpy.zeros((3,),numpy.float) |
| 126 |
counter=0 |
| 127 |
for t in self.getTags(): |
| 128 |
for s in self.getFaultSegments(t)[0]: |
| 129 |
self.__center[:2]+=s[:2] |
| 130 |
counter+=1 |
| 131 |
self.__center/=counter |
| 132 |
return self.__center[:self.getDim()] |
| 133 |
|
| 134 |
def getOrientationOnSurface(self): |
| 135 |
""" |
| 136 |
returns the orientation of the fault system in RAD on the surface around the fault system center |
| 137 |
:rtype: ``float`` |
| 138 |
""" |
| 139 |
if self.__orientation == None: |
| 140 |
center=self.getCenterOnSurface() |
| 141 |
covariant=numpy.zeros((2,2)) |
| 142 |
for t in self.getTags(): |
| 143 |
for s in self.getFaultSegments(t)[0]: |
| 144 |
covariant[0,0]+=(center[0]-s[0])**2 |
| 145 |
covariant[0,1]+=(center[1]-s[1])*(center[0]-s[0]) |
| 146 |
covariant[1,1]+=(center[1]-s[1])**2 |
| 147 |
covariant[1,0]+=(center[1]-s[1])*(center[0]-s[0]) |
| 148 |
e, V=numpy.linalg.eigh(covariant) |
| 149 |
if e[0]>e[1]: |
| 150 |
d=V[:,1] |
| 151 |
else: |
| 152 |
d=V[:,0] |
| 153 |
if abs(d[0])>0.: |
| 154 |
self.__orientation=atan(d[1]/d[0]) |
| 155 |
else: |
| 156 |
self.__orientation=math.pi/2 |
| 157 |
return self.__orientation |
| 158 |
def transform(self, rot=0, shift=numpy.zeros((3,))): |
| 159 |
""" |
| 160 |
applies a shift and a consecutive rotation in the x0x1 plane. |
| 161 |
|
| 162 |
:param rot: rotation angle in RAD |
| 163 |
:type rot: ``float`` |
| 164 |
:param shift: shift vector to be applied before rotation |
| 165 |
:type shift: `numpy.ndarray` of size 2 or 3 |
| 166 |
""" |
| 167 |
if self.getDim() == 2: |
| 168 |
mat=numpy.array([[cos(rot), -sin(rot)], [sin(rot), cos(rot)] ]) |
| 169 |
else: |
| 170 |
mat=numpy.array([[cos(rot), -sin(rot),0.], [sin(rot), cos(rot),0.], [0.,0.,1.] ]) |
| 171 |
for t in self.getTags(): |
| 172 |
top=[] |
| 173 |
for s in self.getFaultSegments(t)[0]: top.append(numpy.dot(mat,(s+shift[:self.getDim()]))) |
| 174 |
if self.getDim() == 3: |
| 175 |
bottom=[] |
| 176 |
for s in self.getFaultSegments(t)[1]: bottom.append(numpy.dot(mat,(s+shift[:self.getDim()]))) |
| 177 |
else: |
| 178 |
bottom=None |
| 179 |
self.addFault(top=top, tag=t, bottom=bottom, w0_offsets=self.getW0Offsets(t), w1_max=-self.getW1Range(t)[0]) |
| 180 |
|
| 181 |
def addFault(self, top, tag=None, bottom=None, w0_offsets=None, w1_max=None): |
| 182 |
""" |
| 183 |
adds a new fault to the fault system. The fault is named by ``tag``. |
| 184 |
|
| 185 |
a fault is subdivided into segments where segment ``i`` is represented by the four corner points |
| 186 |
``v0``=``top[i]``, ``v1``=``bottom[i]``, ``v2``=``bottom[i+1]`` and ``v3``==``top[i+1]``. The segment is parametrized |
| 187 |
by ``w0`` and ``w1`` with ``w0_offsets[i]<=w0<=w0_offsets[i+1]`` and ``-w1_max<=w1<=0``. Moreover |
| 188 |
|
| 189 |
- ``(w0,w1)=(w0_offsets[i] , 0)->v0`` |
| 190 |
- ``(w0,w1)=(w0_offsets[i] , -w1_max)->v1`` |
| 191 |
- ``(w0,w1)=(w0_offsets[i+1], -w1_max)->v2`` |
| 192 |
- ``(w0,w1)=(w0_offsets[i+1], 0)->v3`` |
| 193 |
|
| 194 |
If no ``w0_offsets`` is given, |
| 195 |
|
| 196 |
- ``w0_offsets[0]=0`` |
| 197 |
- ``w0_offsets[i]=w0_offsets[i-1]+length(v0-v4)`` |
| 198 |
|
| 199 |
If no ``w1_max`` is given, the average fault depth (=length(v0-v1)) is used. |
| 200 |
|
| 201 |
:param top: list of coordinates. top[i] and top[i+1] are the top corners of the i-th segment of the fault |
| 202 |
:type top: ``list`` of objects which can be converted to `numpy.ndarray` vector of length 3 (or 2 in the 2D case) |
| 203 |
:param tag: the tag of the fault. If fault ``tag`` already exists it is overwritten. |
| 204 |
:type tag: ``float`` or ``str`` |
| 205 |
:param bottom: list of coordinates. bottom[i] and bottom[i+1] are the bottom corners of the i-th segment of the fault |
| 206 |
:type bottom: ``list`` of objects which can be converted to `numpy.ndarray` vector of length 3 (or 2 in the 2D case) |
| 207 |
:param w0_offsets: ``w0_offsets[i]`` defines the offset of the segment ``i`` in the fault to be used in the parametrization of the fault. If not present the cumulative length of the fault segments is used. |
| 208 |
:type w0_offsets: ``list`` of ``float`` or ``None`` |
| 209 |
:param w1_max: the maximum value used for parametrization of the fault in the depth direction. If not present the mean depth of the fault segments is used. |
| 210 |
:type w1_max: ``float`` |
| 211 |
:note: In the three dimensional case the lists ``bottom`` and ``top`` must have the same length. |
| 212 |
""" |
| 213 |
if tag==None: |
| 214 |
tag=self.NOTAG |
| 215 |
else: |
| 216 |
if self.NOTAG in self.getTags(): |
| 217 |
raise ValueError,'Attempt to add a fault with no tag to a set of existing faults' |
| 218 |
if len(top)<2: |
| 219 |
raise Value,"at least two vertices must be given." |
| 220 |
if self.getDim()==2 and not bottom==None: |
| 221 |
raise ValueError,'Spatial dimension two does not support bottom polygon for faults.' |
| 222 |
if not bottom == None: |
| 223 |
if len(top) != len(bottom): |
| 224 |
raise ValueError,'length of top and bottom polygon must be identical.' |
| 225 |
if w0_offsets != None: |
| 226 |
if len(w0_offsets) != len(top): |
| 227 |
raise ValueError,'expected length of w0_offsets is %s'%(len(top)) |
| 228 |
self.__center=None |
| 229 |
self.__orientation = None |
| 230 |
# |
| 231 |
# translate to numpy |
| 232 |
# |
| 233 |
if self.getDim()==2: |
| 234 |
self.__faults[tag]=([ numpy.array(t[:2],numpy.double) for t in top ],) |
| 235 |
else: |
| 236 |
self.__faults[tag]=( [ numpy.array(t[:3],numpy.double) for t in top ], [ numpy.array(b[:3],numpy.double) for b in bottom ]) |
| 237 |
faults=self.__faults[tag] |
| 238 |
len_faults=len(faults[0]) |
| 239 |
# |
| 240 |
# calculate the depth |
| 241 |
# |
| 242 |
if self.getDim()==2: |
| 243 |
self.__depth[tag]=0. |
| 244 |
else: |
| 245 |
self.__depth[tag]=float(sum([numpy.linalg.norm(faults[1][i]-faults[0][i]) for i in xrange(len_faults)])/len_faults) |
| 246 |
if w1_max==None or self.getDim()==2: w1_max=self.__depth[tag] |
| 247 |
self.__length[tag]=float(sum([numpy.linalg.norm(faults[0][i]-faults[0][i-1]) for i in xrange(1,len_faults)])) |
| 248 |
self.__w1_max[tag]=w1_max |
| 249 |
# |
| 250 |
# |
| 251 |
# |
| 252 |
if w0_offsets!=None: |
| 253 |
self.__offsets[tag]=w0_offsets |
| 254 |
else: |
| 255 |
self.__offsets[tag]=[0.] |
| 256 |
for i in xrange(1,len_faults): |
| 257 |
self.__offsets[tag].append(self.__offsets[tag][-1]+float(numpy.linalg.norm(faults[0][i]-faults[0][i-1]))) |
| 258 |
w0_max=max(self.__offsets[tag]) |
| 259 |
self.__w0_max[tag]=w0_max |
| 260 |
|
| 261 |
def getMaxValue(self,f, tol=sqrt(EPSILON)): |
| 262 |
""" |
| 263 |
returns the maximum value of ``f``, the fault and the location on the fault in fault coordinates. |
| 264 |
|
| 265 |
:param f: a distribution of values |
| 266 |
:type f: `escript.Data` |
| 267 |
:param tol: relative tolerance used to decide if point is on fault |
| 268 |
:type f: ``tol`` |
| 269 |
:return: the maximum value of across all faults in the fault system, the fault tag the maximum is taken, and the coordinates of the location in the coordinates of the fault. The returned fault tag is ``None`` if no points on the fault are available |
| 270 |
""" |
| 271 |
ref=-Lsup(f)*2 |
| 272 |
f_max=ref |
| 273 |
t_max=None |
| 274 |
p_max=None |
| 275 |
x=f.getFunctionSpace().getX() |
| 276 |
for t in self.getTags(): |
| 277 |
p,m=self.getParametrization(x,tag=t, tol=tol) |
| 278 |
loc=((m*f)+(1.-m)*ref).maxGlobalDataPoint() |
| 279 |
f_t=f.getTupleForGlobalDataPoint(*loc)[0] |
| 280 |
if f_t>f_max: |
| 281 |
f_max=f_t |
| 282 |
t_max=t |
| 283 |
p_max=p.getTupleForGlobalDataPoint(*loc)[0] |
| 284 |
|
| 285 |
return f_max, t_max, p_max |
| 286 |
|
| 287 |
def getParametrization(self,x,tag=None, tol=sqrt(EPSILON), outsider=None): |
| 288 |
""" |
| 289 |
returns the parametrization of the fault ``tag`` in the fault system. In fact the values of the parametrization for at given coordinates ``x`` is returned. In addition to the value of the parametrization a mask is returned indicating if the given location is on the fault with given tolerance ``tol``. |
| 290 |
|
| 291 |
Typical usage of the this method is |
| 292 |
|
| 293 |
dom=Domain(..) |
| 294 |
x=dom.getX() |
| 295 |
fs=FaultSystem() |
| 296 |
fs.addFault(tag=3,...) |
| 297 |
p, m=fs.getParametrization(x, outsider=0,tag=3) |
| 298 |
saveDataCSV('x.csv',p=p, x=x, mask=m) |
| 299 |
|
| 300 |
to create a file with the coordinates of the points in ``x`` which are on the fault (as ``mask=m``) together with their location ``p`` in the fault coordinate system. |
| 301 |
|
| 302 |
:param x: location(s) |
| 303 |
:type x: `escript.Data` object or `numpy.ndarray` |
| 304 |
:param tag: the tage of the fault |
| 305 |
:param tol: relative tolerance to check if location is on fault. |
| 306 |
:type tol: ``float`` |
| 307 |
:param outsider: value used for parametrization values outside the fault. If not present an appropriate value is choosen. |
| 308 |
:type outsider: ``float`` |
| 309 |
:return: the coordinates ``x`` in the coordinate system of the fault and a mask indicating coordinates in the fault by 1 (0 elsewhere) |
| 310 |
:rtype: `escript.Data` object or `numpy.ndarray` |
| 311 |
""" |
| 312 |
offsets=self.getW0Offsets(tag) |
| 313 |
w0_range=self.getW0Range(tag)[1]-self.getW0Range(tag)[0] |
| 314 |
if outsider == None: |
| 315 |
outsider=min(self.getW0Range(tag)[0],self.getW0Range(tag)[1])-abs(w0_range)/sqrt(EPSILON) |
| 316 |
|
| 317 |
p=x[0]*0 + outsider |
| 318 |
updated=x[0]*0 |
| 319 |
|
| 320 |
if self.getDim()==2: |
| 321 |
# |
| 322 |
# |
| 323 |
s=self.getFaultSegments(tag)[0] |
| 324 |
for i in xrange(1,len(s)): |
| 325 |
d=s[i]-s[i-1] |
| 326 |
h=x-s[i-1] |
| 327 |
h_l=length(h) |
| 328 |
d_l=length(d) |
| 329 |
if not d_l>0: |
| 330 |
raise ValueError,"segement %s in fault %s has zero length."%(i,tag) |
| 331 |
t=inner(h,d)/d_l**2 |
| 332 |
t=t*whereNonPositive(t-1.)*whereNonNegative(t) |
| 333 |
r=length(h-t*d)/maximum(h_l,d_l) |
| 334 |
m=whereNonPositive(r-tol)*(1.-updated) |
| 335 |
p=(1.-m)*p+m*(offsets[i-1]+(offsets[i]-offsets[i-1])*t) |
| 336 |
updated=wherePositive(updated+m) |
| 337 |
else: |
| 338 |
raise ValueError,"3D is not supported yet." |
| 339 |
|
| 340 |
return p, updated |
| 341 |
|
| 342 |
def __getSignedDistance(self,x,tag=None): |
| 343 |
""" |
| 344 |
returns the signed distance at ``x`` from the fault ``tag`` where the first and the last segments are extended to infinity. |
| 345 |
|
| 346 |
:param x: location(s) |
| 347 |
:type x: `escript.Data` object or `numpy.ndarray` |
| 348 |
:param tag: the tage of the fault |
| 349 |
""" |
| 350 |
p=x[0]*0 |
| 351 |
if self.getDim()==2: |
| 352 |
mat=numpy.array([[0., -1.], [1., 0.] ]) |
| 353 |
# |
| 354 |
# |
| 355 |
s=self.getFaultSegments(tag)[0] |
| 356 |
for i in xrange(1,len(s)): |
| 357 |
d=(s[i]-s[i-1]) |
| 358 |
h=x-s[i-1] |
| 359 |
d_l=length(d) |
| 360 |
if not d_l>0: |
| 361 |
raise ValueError,"segement %s in fault %s has zero length."%(i,tag) |
| 362 |
dt=numpy.dot(mat,d)/d_l |
| 363 |
a=numpy.dot(h,dt) |
| 364 |
p=maximum(a,p) |
| 365 |
else: |
| 366 |
raise ValueError,"3D is not supported yet." |
| 367 |
return p |
| 368 |
|
| 369 |
|
| 370 |
def patchMap(v,v1,v2,v3,v4,x_offset,x_length,depth): |
| 371 |
""" |
| 372 |
maps a patch with corners v1,v2,v3,v4 onto |
| 373 |
the rectangle [x_offset,0] x [x_offset+x_length,-depth] |
| 374 |
""" |
| 375 |
# |
| 376 |
# inspect the mapping v-> w=(w0,w1) with 0<=s<=1 0<=t<=1 |
| 377 |
# then we apply x_offset+s*x_length, -depth*t to map to the final domain |
| 378 |
# |
| 379 |
# we need to find s and t such that |
| 380 |
# |
| 381 |
# v=v1*(1-t)*(1-s)+v2*t*(1-s)+v3*t*s+v4*(1-t)*s |
| 382 |
# =v1*(1-t-s+t*s)+v2*(t-t*s)+v3*t*s+v4*(s-t*s) |
| 383 |
# =v1+t*(v2-v1)+s*(v4-v1)+s*t*(v1-v2+v3-v4) |
| 384 |
# |
| 385 |
# with d4=(v2-v1)-dot(v2-v1,v4-v1)/dot(v4-v1,v4-v1)*(v4-v1) => dot(d4,v2-v1)=0 |
| 386 |
# with d2=(v4-v1)-dot(v2-v1,v4-v1)/dot(v2-v1,v2-v1)*(v2-v1) => dot(d2,v4-v1)=0 |
| 387 |
# |
| 388 |
# dot(v-v1,d4) = s*dot(v4-v1,d4)+s*t*dot(v3-v4,d4) (a) |
| 389 |
# dot(v-v1,d2) = t*dot(v2-v1,d2)+s*t*dot(v3-v2,d2) (b) |
| 390 |
# |
| 391 |
# in addition we set n=cross(d2,d4) => dot(v2-v1,n)=dot(v4-v1,n)=0 |
| 392 |
# |
| 393 |
# dot(v-v1,n)=s*t*dot(v3-v1,n) |
| 394 |
# |
| 395 |
# => s*t=dot(v-v1,n)/dot(v3-v1,n) if dot(v3-v1,n)!=0 |
| 396 |
# otherwise on can set s*t=0 to solve (a) and (b) |
| 397 |
# |
| 398 |
h=v-v1 |
| 399 |
h2=(v2-v1) |
| 400 |
h3=(v3-v1) |
| 401 |
h4=(v4-v1) |
| 402 |
d4=h4-numpy.dot(h2,h4)/numpy.dot(h2,h2)*h2 |
| 403 |
d2=h2-numpy.dot(h2,h4)/numpy.dot(h4,h4)*h4 |
| 404 |
n=numpy.cross(d4,d2) |
| 405 |
h31n=numpy.dot(n,v3-v1) |
| 406 |
if abs(h31n)>0: |
| 407 |
st=n/h31n |
| 408 |
else: |
| 409 |
st=0*n |
| 410 |
s=numpy.dot(h,d4-st*numpy.dot(v3-v4,d4))/numpy.dot(v4-v1,d4) |
| 411 |
t=numpy.dot(h,d2-st*numpy.dot(v3-v2,d2))/numpy.dot(v2-v1,d2) |
| 412 |
return x_offset+s*x_length, -depth*t |
| 413 |
|
Parent Directory
| 
Revision Log