# Contents of /trunk/escript/py_src/faultsystems.py

Revision 2923 - (show annotations)
Thu Feb 4 04:05:36 2010 UTC (9 years, 6 months ago) by jfenwick
File MIME type: text/x-python
File size: 24968 byte(s)
```Bringing non-release specific things from stage3.1 r2922 back to trunk

```
 1 ######################################################## 2 # 3 # Copyright (c) 2003-2010 by University of Queensland 4 # Earth Systems Science Computational Center (ESSCC) 5 6 # 7 # Primary Business: Queensland, Australia 8 # Licensed under the Open Software License version 3.0 9 10 # 11 ######################################################## 12 13 __copyright__="""Copyright (c) 2003-2010 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 19 __url__= 20 21 from esys.escript import * 22 from esys.escript.pdetools import Locator 23 import numpy 24 import math 25 26 class FaultSystem: 27 """ 28 The FaultSystem class defines a system of faults in the Earth's crust. 29 30 A fault system is defined by set of faults index by a tag. Each fault is defined by a starting point V0 and a list of 31 strikes ``strike`` and length ``l``. The strikes and the length is used to define a polyline with points ``V[i]`` such that 32 33 - ``V[0]=V0`` 34 - ``V[i]=V[i]+ls[i]*array(cos(strike[i]),sin(strike[i]),0)'' 35 36 So ``strike`` defines the angle between the direction of the fault segment and the x0 axis. ''ls[i]``==0 is allowed. 37 38 In case of a 3D model a fault plane is defined through a dip and depth. 39 40 The class provides a mechanism to parametrise each fault with the domain [0,w0_max] x [0, w1_max] (to [0,w0_max] in the 2D case). 41 """ 42 NOTAG="__NOTAG__" 43 MIN_DEPTH_ANGLE=0.1 44 def __init__(self,dim=3): 45 """ 46 Sets up the fault system 47 48 :param dim: spatial dimension 49 :type dim: ``int`` of value 2 or 3 50 """ 51 if not (dim == 2 or dim == 3): 52 raise ValueError,"only dimension2 2 and 3 are supported." 53 self.__dim=dim 54 self.__top={} 55 self.__ls={} 56 self.__strikes={} 57 self.__strike_vectors={} 58 self.__medDepth={} 59 self.__total_length={} 60 if dim ==2: 61 self.__depths=None 62 self.__depth_vectors=None 63 self.__dips=None 64 self.__bottom=None 65 self.__normals=None 66 else: 67 self.__depths={} 68 self.__depth_vectors={} 69 self.__dips={} 70 self.__bottom={} 71 self.__normals={} 72 self.__offsets={} 73 self.__w1_max={} 74 self.__w0_max={} 75 self.__center=None 76 self.__orientation = None 77 def getStart(self,tag=None): 78 """ 79 returns the starting point of fault ``tag`` 80 :rtype: `numpy.ndarray'. 81 """ 82 return self.getTopPolyline(tag)[0] 83 84 def getTags(self): 85 """ 86 returns a list of the tags used by the fault system 87 :rtype: ``list`` 88 """ 89 return self.__top.keys() 90 def getDim(self): 91 """ 92 returns the spatial dimension 93 :rtype: ``int`` 94 """ 95 return self.__dim 96 97 def getTopPolyline(self, tag=None): 98 """ 99 returns the polyline used to describe fault tagged by ``tag`` 100 101 :param tag: the tag of the fault 102 :type tag: ``float`` or ``str`` 103 :return: the list of vertices defining the top of the fault. The coordinates are `numpy.ndarray'. 104 """ 105 if tag==None: tag=self.NOTAG 106 return self.__top[tag] 107 def getStrikes(self, tag=None): 108 """ 109 returns the strike of the segements in fault ``tag`` 110 :rtype: ``list`` of ``float`` 111 """ 112 if tag==None: tag=self.NOTAG 113 return self.__strikes[tag] 114 def getStrikeVectors(self, tag=None): 115 """ 116 returns the strike vectors of fault ``tag`` 117 :rtype: ``list`` of `numpy.ndarray'. 118 """ 119 if tag==None: tag=self.NOTAG 120 return self.__strike_vectors[tag] 121 def getLengths(self, tag=None): 122 """ 123 returns the lengths of segments in fault ``tag`` 124 :rtype: ``list`` of `float`` 125 """ 126 if tag==None: tag=self.NOTAG 127 return self.__ls[tag] 128 129 def getTotalLength(self, tag=None): 130 """ 131 returns the total unrolled length of fault ``tag`` 132 :rtype: `float`` 133 """ 134 if tag==None: tag=self.NOTAG 135 return self.__total_length[tag] 136 137 def getMediumDepth(self,tag=None): 138 """ 139 returns the medium depth of fault ``tag`` 140 :rtype: ``float`` 141 """ 142 if tag==None: tag=self.NOTAG 143 return self.__medDepth[tag] 144 145 def getDips(self, tag=None): 146 """ 147 returns the list of the dips of the segements in fault ``tag`` 148 :param tag: the tag of the fault 149 :type tag: ``float`` or ``str`` 150 :return: the list of segment dips. In the 2D case None is returned. 151 """ 152 if tag==None: tag=self.NOTAG 153 if self.getDim()==3: 154 return self.__dips[tag] 155 else: 156 return None 157 158 def getBottomPolyline(self, tag=None): 159 """ 160 returns the list of the vertices defining the bottom of the fault ``tag`` 161 :param tag: the tag of the fault 162 :type tag: ``float`` or ``str`` 163 :return: the list of vertices. In the 2D case None is returned. 164 """ 165 if tag==None: tag=self.NOTAG 166 if self.getDim()==3: 167 return self.__bottom[tag] 168 else: 169 return None 170 171 def getSegmentNormals(self, tag=None): 172 """ 173 returns the list of the normals of the segments in fault ``tag`` 174 :param tag: the tag of the fault 175 :type tag: ``float`` or ``str`` 176 :return: the list of vectors normal to the segments. In the 2D case None is returned. 177 """ 178 if tag==None: tag=self.NOTAG 179 if self.getDim()==3: 180 return self.__normals[tag] 181 else: 182 return None 183 184 def getDepthVectors(self, tag=None): 185 """ 186 returns the list of the depth vector at top vertices in fault ``tag``. 187 :param tag: the tag of the fault 188 :type tag: ``float`` or ``str`` 189 :return: the list of segment depths. In the 2D case None is returned. 190 """ 191 if tag==None: tag=self.NOTAG 192 if self.getDim()==3: 193 return self.__depth_vectors[tag] 194 else: 195 return None 196 def getDepths(self, tag=None): 197 """ 198 returns the list of the depths of the segements in fault ``tag``. 199 :param tag: the tag of the fault 200 :type tag: ``float`` or ``str`` 201 :return: the list of segment depths. In the 2D case None is returned. 202 """ 203 if tag==None: tag=self.NOTAG 204 if self.getDim()==3: 205 return self.__depths[tag] 206 else: 207 return None 208 209 def getW0Range(self,tag=None): 210 """ 211 returns the range of the parameterization in ``w0`` 212 :rtype: two ``float`` 213 """ 214 return self.getW0Offsets(tag)[0], self.getW0Offsets(tag)[-1] 215 216 def getW1Range(self,tag=None): 217 """ 218 returns the range of the parameterization in ``w1`` 219 :rtype: two ``float`` 220 """ 221 if tag==None: tag=self.NOTAG 222 return -self.__w1_max[tag],0 223 224 def getW0Offsets(self, tag=None): 225 """ 226 returns the offsets for the parametrization of fault ``tag``. 227 228 :return: the offsets in the parametrization 229 :rtype: ``list`` of ``float`` 230 """ 231 if tag==None: tag=self.NOTAG 232 return self.__offsets[tag] 233 234 235 def getCenterOnSurface(self): 236 """ 237 returns the center point of the fault system at the surface 238 :rtype: `numpy.ndarray` 239 """ 240 if self.__center == None: 241 self.__center=numpy.zeros((3,),numpy.float) 242 counter=0 243 for t in self.getTags(): 244 for s in self.getTopPolyline(t): 245 self.__center[:2]+=s[:2] 246 counter+=1 247 self.__center/=counter 248 return self.__center[:self.getDim()] 249 250 def getOrientationOnSurface(self): 251 """ 252 returns the orientation of the fault system in RAD on the surface around the fault system center 253 :rtype: ``float`` 254 """ 255 if self.__orientation == None: 256 center=self.getCenterOnSurface() 257 covariant=numpy.zeros((2,2)) 258 for t in self.getTags(): 259 for s in self.getTopPolyline(t): 260 covariant[0,0]+=(center[0]-s[0])**2 261 covariant[0,1]+=(center[1]-s[1])*(center[0]-s[0]) 262 covariant[1,1]+=(center[1]-s[1])**2 263 covariant[1,0]+=(center[1]-s[1])*(center[0]-s[0]) 264 e, V=numpy.linalg.eigh(covariant) 265 if e[0]>e[1]: 266 d=V[:,0] 267 else: 268 d=V[:,1] 269 if abs(d[0])>0.: 270 self.__orientation=atan(d[1]/d[0]) 271 else: 272 self.__orientation=math.pi/2 273 return self.__orientation 274 def transform(self, rot=0, shift=numpy.zeros((3,))): 275 """ 276 applies a shift and a consecutive rotation in the x0x1 plane. 277 278 :param rot: rotation angle in RAD 279 :type rot: ``float`` 280 :param shift: shift vector to be applied before rotation 281 :type shift: `numpy.ndarray` of size 2 or 3 282 """ 283 if self.getDim() == 2: 284 mat=numpy.array([[cos(rot), -sin(rot)], [sin(rot), cos(rot)] ]) 285 else: 286 mat=numpy.array([[cos(rot), -sin(rot),0.], [sin(rot), cos(rot),0.], [0.,0.,1.] ]) 287 288 for t in self.getTags(): 289 strikes=[ s+ rot for s in self.getStrikes(t) ] 290 V0=self.getStart(t) 291 292 self.addFault(strikes = [ s+ rot for s in self.getStrikes(t) ], \ 293 ls = self.getLengths(t), \ 294 V0=numpy.dot(mat,self.getStart(t)+shift), \ 295 tag =t, \ 296 dips=self.getDips(t),\ 297 depths=self.getDepths(t), \ 298 w0_offsets=self.getW0Offsets(t), \ 299 w1_max=-self.getW1Range(t)[0]) 300 301 def addFault(self, strikes, ls, V0=[0.,0.,0.],tag=None, dips=None, depths= None, w0_offsets=None, w1_max=None): 302 """ 303 adds a new fault to the fault system. The fault is named by ``tag``. 304 305 The fault is defined by a starting point V0 and a list of ``strikes`` and length ``ls``. The strikes and the length 306 is used to define a polyline with points ``V[i]`` such that 307 308 - ``V[0]=V0`` 309 - ``V[i]=V[i]+ls[i]*array(cos(strikes[i]),sin(strikes[i]),0)'' 310 311 So ``strikes`` defines the angle between the direction of the fault segment and the x0 axis. In 3D ''ls[i]``==0 is allowed. 312 313 In case of a 3D model a fault plane is defined through a dip ``dips`` and depth ``depths``. 314 From the dip and the depth the polyline ``bottom`` of the bottom of the fault is computed. 315 316 317 Each segment in the fault is decribed by the for vertices ``v0=top[i]``, ``v1``==``top[i+1]``, ``v2=bottom[i]`` and ``v3=bottom[i+1]`` 318 The segment is parametrized by ``w0`` and ``w1`` with ``w0_offsets[i]<=w0<=w0_offsets[i+1]`` and ``-w1_max<=w1<=0``. Moreover 319 320 - ``(w0,w1)=(w0_offsets[i] , 0)->v0`` 321 - ``(w0,w1)=(w0_offsets[i+1], 0)->v1`` 322 - ``(w0,w1)=(w0_offsets[i] , -w1_max)->v2`` 323 - ``(w0,w1)=(w0_offsets[i+1], -w1_max)->v3`` 324 325 If no ``w0_offsets`` is given, 326 327 - ``w0_offsets[0]=0`` 328 - ``w0_offsets[i]=w0_offsets[i-1]+L[i]`` 329 330 where ``L[i]`` is the length of the segments on the top in 2D and in the middle of the segment in 3D. 331 332 If no ``w1_max`` is given, the average fault depth is used. 333 334 335 :param strikes: list of strikes. This is the angle of the fault segment direction with x0 axis. Right hand rule applies. 336 :type strikes: ``list`` of ``float`` 337 :param ls: list of fault lengths. In the case of a 3D fault a segment may have length 0. 338 :type ls: ``list`` of ``float`` 339 :param V0: start point of the fault 340 :type V0: ``list`` or `numpy.ndarray' with 2 or 3 components. ``V0[2]`` must be zero. 341 :param tag: the tag of the fault. If fault ``tag`` already exists it is overwritten. 342 :type tag: ``float`` or ``str`` 343 :param dips: list of dip angles. Right hand rule around strike direction applies. 344 :type dips: ``list`` of ``float`` 345 :param depths: list of segment depth. Value mut be positive in the 3D case. 346 :type depths: ``list`` of ``float`` 347 :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. 348 :type w0_offsets: ``list`` of ``float`` or ``None`` 349 :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. 350 :type w1_max: ``float`` 351 :note: In the three dimensional case the lists ``dip`` and ``top`` must have the same length. 352 """ 353 if tag==None: 354 tag=self.NOTAG 355 else: 356 if self.NOTAG in self.getTags(): 357 raise ValueError,'Attempt to add a fault with no tag to a set of existing faults' 358 if not isinstance(strikes, list): strikes=[strikes, ] 359 n_segs=len(strikes) 360 if not isinstance(ls, list): ls=[ ls for i in xrange(n_segs) ] 361 if not n_segs==len(ls): 362 raise ValueError,"number of strike direction and length must match." 363 if len(V0)>2: 364 if abs(V0[2])>0: raise Value,"start point needs to be surface (3rd component ==0)" 365 if self.getDim()==2 and not (dips==None and depths == None) : 366 raise ValueError,'Spatial dimension two does not support dip and depth for faults.' 367 if not dips == None: 368 if not isinstance(dips, list): dips=[dips for i in xrange(n_segs) ] 369 if n_segs != len(dips): 370 raise ValueError,'length of dips must be one less then the length of top.' 371 if not depths == None: 372 if not isinstance(depths, list): depths=[depths for i in xrange(n_segs+1) ] 373 if n_segs+1 != len(depths): 374 raise ValueError,'length of depths must be one less then the length of top.' 375 if w0_offsets != None: 376 if len(w0_offsets) != n_segs+1: 377 raise ValueError,'expected length of w0_offsets is %s'%(n_segs) 378 self.__center=None 379 self.__orientation = None 380 # 381 # in the 2D case we don't allow zero length: 382 # 383 if self.getDim() == 2: 384 for l in ls: 385 if l<=0: raise ValueError,"length must be positive" 386 else: 387 for l in ls: 388 if l<0: raise ValueError,"length must be non-negative" 389 for i in xrange(n_segs+1): 390 if depths[i]<0: raise ValueError,"negative depth." 391 # 392 # translate start point to numarray 393 # 394 V0= numpy.array(V0[:self.getDim()],numpy.double) 395 # 396 # set strike vectors: 397 # 398 strike_vectors=[] 399 top_polyline=[V0] 400 total_length=0 401 for i in xrange(n_segs): 402 v=numpy.zeros((self.getDim(),)) 403 v[0]=cos(strikes[i]) 404 v[1]=sin(strikes[i]) 405 strike_vectors.append(v) 406 top_polyline.append(top_polyline[-1]+ls[i]*v) 407 total_length+=ls[i] 408 # 409 # normal and depth direction 410 # 411 if self.getDim()==3: 412 normals=[] 413 for i in xrange(n_segs): 414 normals.append(numpy.array([sin(dips[i])*strike_vectors[i][1],-sin(dips[i])*strike_vectors[i][0], cos(dips[i])]) ) 415 416 d=numpy.cross(strike_vectors[0],normals[0]) 417 if d[2]>0: 418 f=-1 419 else: 420 f=1 421 depth_vectors=[f*depths[0]*d/numpy.linalg.norm(d) ] 422 for i in xrange(1,n_segs): 423 d=-numpy.cross(normals[i-1],normals[i]) 424 d_l=numpy.linalg.norm(d) 425 if d_l<=0: 426 d=numpy.cross(strike_vectors[i],normals[i]) 427 d_l=numpy.linalg.norm(d) 428 else: 429 for L in [ strike_vectors[i], strike_vectors[i-1]]: 430 if numpy.linalg.norm(numpy.cross(L,d)) <= self.MIN_DEPTH_ANGLE * numpy.linalg.norm(L) * d_l: 431 raise ValueError,"%s-th depth vector %s too flat."%(i, d) 432 if d[2]>0: 433 f=-1 434 else: 435 f=1 436 depth_vectors.append(f*d*depths[i]/d_l) 437 d=numpy.cross(strike_vectors[n_segs-1],normals[n_segs-1]) 438 if d[2]>0: 439 f=-1 440 else: 441 f=1 442 depth_vectors.append(f*depths[n_segs]*d/numpy.linalg.norm(d)) 443 bottom_polyline=[ top_polyline[i]+depth_vectors[i] for i in xrange(n_segs+1) ] 444 # 445 # calculate offsets if required: 446 # 447 if w0_offsets==None: 448 w0_offsets=[0.] 449 for i in xrange(n_segs): 450 if self.getDim()==3: 451 w0_offsets.append(w0_offsets[-1]+(float(numpy.linalg.norm(bottom_polyline[i+1]-bottom_polyline[i]))+ls[i])/2.) 452 else: 453 w0_offsets.append(w0_offsets[-1]+ls[i]) 454 w0_max=max(w0_offsets) 455 if self.getDim()==3: 456 self.__normals[tag]=normals 457 self.__depth_vectors[tag]=depth_vectors 458 self.__depths[tag]=depths 459 self.__dips[tag]=dips 460 self.__bottom[tag]=bottom_polyline 461 self.__ls[tag]=ls 462 self.__strikes[tag]=strikes 463 self.__strike_vectors[tag]=strike_vectors 464 self.__top[tag]=top_polyline 465 self.__total_length[tag]=total_length 466 self.__offsets[tag]=w0_offsets 467 468 if self.getDim()==2: 469 self.__medDepth[tag]=0. 470 else: 471 self.__medDepth[tag]=sum([ numpy.linalg.norm(v) for v in depth_vectors])/len(depth_vectors) 472 if w1_max==None or self.getDim()==2: w1_max=self.__medDepth[tag] 473 self.__w0_max[tag]=w0_max 474 self.__w1_max[tag]=w1_max 475 476 def getMaxValue(self,f, tol=sqrt(EPSILON)): 477 """ 478 returns the tag of the fault of where ``f`` takes the maximum value and a `Locator` object which can be used to collect values from `Data` class objects at the location where the minimum is taken. 479 480 :param f: a distribution of values 481 :type f: `escript.Data` 482 :param tol: relative tolerance used to decide if point is on fault 483 :type f: ``tol`` 484 :return: the fault tag the maximum is taken, and a `Locator` object to collect the value at location of maximum value. 485 """ 486 ref=-Lsup(f)*2 487 f_max=ref 488 t_max=None 489 loc_max=None 490 x=f.getFunctionSpace().getX() 491 for t in self.getTags(): 492 p,m=self.getParametrization(x,tag=t, tol=tol) 493 loc=((m*f)+(1.-m)*ref).maxGlobalDataPoint() 494 f_t=f.getTupleForGlobalDataPoint(*loc)[0] 495 if f_t>f_max: 496 f_max=f_t 497 t_max=t 498 loc_max=loc 499 500 if loc_max == None: 501 return None, None 502 else: 503 return t_max, Locator(x.getFunctionSpace(),x.getTupleForGlobalDataPoint(*loc_max)) 504 505 def getMinValue(self,f, tol=sqrt(EPSILON)): 506 """ 507 returns the tag of the fault of where ``f`` takes the minimum value and a `Locator` object which can be used to collect values from `Data` class objects at the location where the minimum is taken. 508 509 :param f: a distribution of values 510 :type f: `escript.Data` 511 :param tol: relative tolerance used to decide if point is on fault 512 :type f: ``tol`` 513 :return: the fault tag the minimum is taken, and a `Locator` object to collect the value at location of minimum value. 514 """ 515 ref=Lsup(f)*2 516 f_min=ref 517 t_min=None 518 loc_min=None 519 x=f.getFunctionSpace().getX() 520 for t in self.getTags(): 521 p,m=self.getParametrization(x,tag=t, tol=tol) 522 loc=((m*f)+(1.-m)*ref).minGlobalDataPoint() 523 f_t=f.getTupleForGlobalDataPoint(*loc)[0] 524 if f_t