/[escript]/trunk/escript/py_src/faultsystems.py
ViewVC logotype

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

Parent Directory Parent Directory | Revision Log Revision Log


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 # 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-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 http://www.opensource.org/licenses/osl-3.0.php"""
19 __url__="https://launchpad.net/escript-finley"
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<f_min:
525 f_min=f_t
526 t_min=t
527 loc_min=loc
528
529 if loc_min == None:
530 return None, None
531 else:
532 return t_min, Locator(x.getFunctionSpace(),x.getTupleForGlobalDataPoint(*loc_min))
533
534 def getParametrization(self,x,tag=None, tol=sqrt(EPSILON), outsider=None):
535 """
536 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``.
537
538 Typical usage of the this method is
539
540 dom=Domain(..)
541 x=dom.getX()
542 fs=FaultSystem()
543 fs.addFault(tag=3,...)
544 p, m=fs.getParametrization(x, outsider=0,tag=3)
545 saveDataCSV('x.csv',p=p, x=x, mask=m)
546
547 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.
548
549 :param x: location(s)
550 :type x: `escript.Data` object or `numpy.ndarray`
551 :param tag: the tag of the fault
552 :param tol: relative tolerance to check if location is on fault.
553 :type tol: ``float``
554 :param outsider: value used for parametrization values outside the fault. If not present an appropriate value is choosen.
555 :type outsider: ``float``
556 :return: the coordinates ``x`` in the coordinate system of the fault and a mask indicating coordinates in the fault by 1 (0 elsewhere)
557 :rtype: `escript.Data` object or `numpy.ndarray`
558 """
559 offsets=self.getW0Offsets(tag)
560 w1_range=self.getW1Range(tag)
561 w0_range=self.getW0Range(tag)[1]-self.getW0Range(tag)[0]
562 if outsider == None:
563 outsider=min(self.getW0Range(tag)[0],self.getW0Range(tag)[1])-abs(w0_range)/sqrt(EPSILON)
564
565 if isinstance(x,list): x=numpy.array(x, numpy.double)
566 updated=x[0]*0
567
568 if self.getDim()==2:
569 #
570 #
571 p=x[0]*0 + outsider
572 top=self.getTopPolyline(tag)
573 for i in xrange(1,len(top)):
574 d=top[i]-top[i-1]
575 h=x-top[i-1]
576 h_l=length(h)
577 d_l=length(d)
578 s=inner(h,d)/d_l**2
579 s=s*whereNonPositive(s-1.-tol)*whereNonNegative(s+tol)
580 m=whereNonPositive(length(h-s*d)-tol*maximum(h_l,d_l))*(1.-updated)
581 p=(1.-m)*p+m*(offsets[i-1]+(offsets[i]-offsets[i-1])*s)
582 updated=wherePositive(updated+m)
583 else:
584 p=x[:2]*0 + outsider
585 top=self.getTopPolyline(tag)
586 bottom=self.getBottomPolyline(tag)
587 n=self.getSegmentNormals(tag)
588 for i in xrange(len(top)-1):
589 h=x-top[i]
590 R=top[i+1]-top[i]
591 r=bottom[i+1]-bottom[i]
592 D0=bottom[i]-top[i]
593 D1=bottom[i+1]-top[i+1]
594 s_upper=matrix_mult(numpy.linalg.pinv(numpy.vstack((R,D1)).T),h)
595 s_lower=matrix_mult(numpy.linalg.pinv(numpy.vstack((r,D0)).T),h)
596 m_ul=wherePositive(s_upper[0]-s_upper[1])
597 s=s_upper*m_ul+s_lower*(1-m_ul)
598 s0=s[0]
599 s1=s[1]
600 m=whereNonNegative(s0+tol)*whereNonPositive(s0-1.-tol)*whereNonNegative(s1+tol)*whereNonPositive(s1-1.-tol)
601 s0=s0*m
602 s1=s1*m
603 atol=tol*maximum(length(h),length(top[i]-bottom[i+1]))
604 m=whereNonPositive(length(h-s0*R-s1*D1)*m_ul+(1-m_ul)*length(h-s0*r-s1*D0)-atol)
605 p[0]=(1.-m)*p[0]+m*(offsets[i]+(offsets[i+1]-offsets[i])*s0)
606 p[1]=(1.-m)*p[1]+m*(w1_range[1]+(w1_range[0]-w1_range[1])*s1)
607 updated=wherePositive(updated+m)
608
609 return p, updated
610
611 def getSideAndDistance(self,x,tag=None):
612 """
613 returns the side and the distance at ``x`` from the fault ``tag``.
614
615 :param x: location(s)
616 :type x: `escript.Data` object or `numpy.ndarray`
617 :param tag: the tag of the fault
618 :return: the side of ``x`` (positive means to the right of the fault, negative to the left) and the distance to the fault. Note that a value zero for the side means that that the side is undefined.
619 """
620 d=None
621 side=None
622 if self.getDim()==2:
623 mat=numpy.array([[0., 1.], [-1., 0.] ])
624 s=self.getTopPolyline(tag)
625 for i in xrange(1,len(s)):
626 q=(s[i]-s[i-1])
627 h=x-s[i-1]
628 q_l=length(q)
629 qt=matrixmult(mat,q) # orthogonal direction
630 t=inner(q,h)/q_l**2
631 t=maximum(minimum(t,1,),0.)
632 p=h-t*q
633 dist=length(p)
634 lside=sign(inner(p,qt))
635 if d == None:
636 d=dist
637 side=lside
638 else:
639 m=whereNegative(d-dist)
640 m2=wherePositive(whereZero(abs(lside))+m)
641 d=dist*(1-m)+d*m
642 side=lside*(1-m2)+side*m2
643 else:
644 ns=self.getSegmentNormals(tag)
645 top=self.getTopPolyline(tag)
646 bottom=self.getBottomPolyline(tag)
647 for i in xrange(len(top)-1):
648 h=x-top[i]
649 R=top[i+1]-top[i]
650 r=bottom[i+1]-bottom[i]
651 D0=bottom[i]-top[i]
652 D1=bottom[i+1]-top[i+1]
653 s_upper=matrix_mult(numpy.linalg.pinv(numpy.vstack((R,D1)).T),h)
654 s_lower=matrix_mult(numpy.linalg.pinv(numpy.vstack((r,D0)).T),h)
655 m_ul=wherePositive(s_upper[0]-s_upper[1])
656 s=s_upper*m_ul+s_lower*(1-m_ul)
657 s=maximum(minimum(s,1.),0)
658 p=h-(m_ul*R+(1-m_ul)*r)*s[0]-(m_ul*D1+(1-m_ul)*D0)*s[1]
659 dist=length(p)
660 lside=sign(inner(p,ns[i]))
661 if d == None:
662 d=dist
663 side=lside
664 else:
665 m=whereNegative(d-dist)
666 m2=wherePositive(whereZero(abs(lside))+m)
667 d=dist*(1-m)+d*m
668 side=lside*(1-m2)+side*m2
669
670 return side, d
671

  ViewVC Help
Powered by ViewVC 1.1.26