/[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 2675 - (show annotations)
Mon Sep 21 02:32:17 2009 UTC (11 years, 1 month ago) by gross
File MIME type: text/x-python
File size: 24801 byte(s)
getMinValue added to the fault system.
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 by a starting point V0 and a list of
30 strikes ``strike`` and length ``l``. The strikes and the length is used to define a polyline with points ``V[i]`` such that
31
32 - ``V[0]=V0``
33 - ``V[i]=V[i]+ls[i]*array(cos(strike[i]),sin(strike[i]),0)''
34
35 So ``strike`` defines the angle between the direction of the fault segment and the x0 axis. ''ls[i]``==0 is allowed.
36
37 In case of a 3D model a fault plane is defined through a dip and depth.
38
39 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).
40 """
41 NOTAG="__NOTAG__"
42 MIN_DEPTH_ANGLE=0.1
43 def __init__(self,dim=3):
44 """
45 Sets up the fault system
46
47 :param dim: spatial dimension
48 :type dim: ``int`` of value 2 or 3
49 """
50 if not (dim == 2 or dim == 3):
51 raise ValueError,"only dimension2 2 and 3 are supported."
52 self.__dim=dim
53 self.__top={}
54 self.__ls={}
55 self.__strikes={}
56 self.__strike_vectors={}
57 self.__medDepth={}
58 self.__total_length={}
59 if dim ==2:
60 self.__depths=None
61 self.__depth_vectors=None
62 self.__dips=None
63 self.__bottom=None
64 self.__normals=None
65 else:
66 self.__depths={}
67 self.__depth_vectors={}
68 self.__dips={}
69 self.__bottom={}
70 self.__normals={}
71 self.__offsets={}
72 self.__w1_max={}
73 self.__w0_max={}
74 self.__center=None
75 self.__orientation = None
76 def getStart(self,tag=None):
77 """
78 returns the starting point of fault ``tag``
79 :rtype: `numpy.ndarray'.
80 """
81 return self.getTopPolyline(tag)[0]
82
83 def getTags(self):
84 """
85 returns a list of the tags used by the fault system
86 :rtype: ``list``
87 """
88 return self.__top.keys()
89 def getDim(self):
90 """
91 returns the spatial dimension
92 :rtype: ``int``
93 """
94 return self.__dim
95
96 def getTopPolyline(self, tag=None):
97 """
98 returns the polyline used to describe fault tagged by ``tag``
99
100 :param tag: the tag of the fault
101 :type tag: ``float`` or ``str``
102 :return: the list of vertices defining the top of the fault. The coordinates are `numpy.ndarray'.
103 """
104 if tag==None: tag=self.NOTAG
105 return self.__top[tag]
106 def getStrikes(self, tag=None):
107 """
108 returns the strike of the segements in fault ``tag``
109 :rtype: ``list`` of ``float``
110 """
111 if tag==None: tag=self.NOTAG
112 return self.__strikes[tag]
113 def getStrikeVectors(self, tag=None):
114 """
115 returns the strike vectors of fault ``tag``
116 :rtype: ``list`` of `numpy.ndarray'.
117 """
118 if tag==None: tag=self.NOTAG
119 return self.__strike_vectors[tag]
120 def getLengths(self, tag=None):
121 """
122 returns the lengths of segments in fault ``tag``
123 :rtype: ``list`` of `float``
124 """
125 if tag==None: tag=self.NOTAG
126 return self.__ls[tag]
127
128 def getTotalLength(self, tag=None):
129 """
130 returns the total unrolled length of fault ``tag``
131 :rtype: `float``
132 """
133 if tag==None: tag=self.NOTAG
134 return self.__total_length[tag]
135
136 def getMediumDepth(self,tag=None):
137 """
138 returns the medium depth of fault ``tag``
139 :rtype: ``float``
140 """
141 if tag==None: tag=self.NOTAG
142 return self.__medDepth[tag]
143
144 def getDips(self, tag=None):
145 """
146 returns the list of the dips of the segements in fault ``tag``
147 :param tag: the tag of the fault
148 :type tag: ``float`` or ``str``
149 :return: the list of segment dips. In the 2D case None is returned.
150 """
151 if tag==None: tag=self.NOTAG
152 if self.getDim()==3:
153 return self.__dips[tag]
154 else:
155 return None
156
157 def getBottomPolyline(self, tag=None):
158 """
159 returns the list of the vertices defining the bottom of the fault ``tag``
160 :param tag: the tag of the fault
161 :type tag: ``float`` or ``str``
162 :return: the list of vertices. In the 2D case None is returned.
163 """
164 if tag==None: tag=self.NOTAG
165 if self.getDim()==3:
166 return self.__bottom[tag]
167 else:
168 return None
169
170 def getSegmentNormals(self, tag=None):
171 """
172 returns the list of the normals of the segments in fault ``tag``
173 :param tag: the tag of the fault
174 :type tag: ``float`` or ``str``
175 :return: the list of vectors normal to the segments. In the 2D case None is returned.
176 """
177 if tag==None: tag=self.NOTAG
178 if self.getDim()==3:
179 return self.__normals[tag]
180 else:
181 return None
182
183 def getDepthVectors(self, tag=None):
184 """
185 returns the list of the depth vector at top vertices in fault ``tag``.
186 :param tag: the tag of the fault
187 :type tag: ``float`` or ``str``
188 :return: the list of segment depths. In the 2D case None is returned.
189 """
190 if tag==None: tag=self.NOTAG
191 if self.getDim()==3:
192 return self.__depth_vectors[tag]
193 else:
194 return None
195 def getDepths(self, tag=None):
196 """
197 returns the list of the depths of the segements in fault ``tag``.
198 :param tag: the tag of the fault
199 :type tag: ``float`` or ``str``
200 :return: the list of segment depths. In the 2D case None is returned.
201 """
202 if tag==None: tag=self.NOTAG
203 if self.getDim()==3:
204 return self.__depths[tag]
205 else:
206 return None
207
208 def getW0Range(self,tag=None):
209 """
210 returns the range of the parameterization in ``w0``
211 :rtype: two ``float``
212 """
213 return self.getW0Offsets(tag)[0], self.getW0Offsets(tag)[-1]
214
215 def getW1Range(self,tag=None):
216 """
217 returns the range of the parameterization in ``w1``
218 :rtype: two ``float``
219 """
220 if tag==None: tag=self.NOTAG
221 return -self.__w1_max[tag],0
222
223 def getW0Offsets(self, tag=None):
224 """
225 returns the offsets for the parametrization of fault ``tag``.
226
227 :return: the offsets in the parametrization
228 :rtype: ``list`` of ``float``
229 """
230 if tag==None: tag=self.NOTAG
231 return self.__offsets[tag]
232
233
234 def getCenterOnSurface(self):
235 """
236 returns the center point of the fault system at the surface
237 :rtype: `numpy.ndarray`
238 """
239 if self.__center == None:
240 self.__center=numpy.zeros((3,),numpy.float)
241 counter=0
242 for t in self.getTags():
243 for s in self.getTopPolyline(t):
244 self.__center[:2]+=s[:2]
245 counter+=1
246 self.__center/=counter
247 return self.__center[:self.getDim()]
248
249 def getOrientationOnSurface(self):
250 """
251 returns the orientation of the fault system in RAD on the surface around the fault system center
252 :rtype: ``float``
253 """
254 if self.__orientation == None:
255 center=self.getCenterOnSurface()
256 covariant=numpy.zeros((2,2))
257 for t in self.getTags():
258 for s in self.getTopPolyline(t):
259 covariant[0,0]+=(center[0]-s[0])**2
260 covariant[0,1]+=(center[1]-s[1])*(center[0]-s[0])
261 covariant[1,1]+=(center[1]-s[1])**2
262 covariant[1,0]+=(center[1]-s[1])*(center[0]-s[0])
263 e, V=numpy.linalg.eigh(covariant)
264 if e[0]>e[1]:
265 d=V[:,0]
266 else:
267 d=V[:,1]
268 if abs(d[0])>0.:
269 self.__orientation=atan(d[1]/d[0])
270 else:
271 self.__orientation=math.pi/2
272 return self.__orientation
273 def transform(self, rot=0, shift=numpy.zeros((3,))):
274 """
275 applies a shift and a consecutive rotation in the x0x1 plane.
276
277 :param rot: rotation angle in RAD
278 :type rot: ``float``
279 :param shift: shift vector to be applied before rotation
280 :type shift: `numpy.ndarray` of size 2 or 3
281 """
282 if self.getDim() == 2:
283 mat=numpy.array([[cos(rot), -sin(rot)], [sin(rot), cos(rot)] ])
284 else:
285 mat=numpy.array([[cos(rot), -sin(rot),0.], [sin(rot), cos(rot),0.], [0.,0.,1.] ])
286
287 for t in self.getTags():
288 strikes=[ s+ rot for s in self.getStrikes(t) ]
289 V0=self.getStart(t)
290
291 self.addFault(strikes = [ s+ rot for s in self.getStrikes(t) ], \
292 ls = self.getLengths(t), \
293 V0=numpy.dot(mat,self.getStart(t)+shift), \
294 tag =t, \
295 dips=self.getDips(t),\
296 depths=self.getDepths(t), \
297 w0_offsets=self.getW0Offsets(t), \
298 w1_max=-self.getW1Range(t)[0])
299
300 def addFault(self, strikes, ls, V0=[0.,0.,0.],tag=None, dips=None, depths= None, w0_offsets=None, w1_max=None):
301 """
302 adds a new fault to the fault system. The fault is named by ``tag``.
303
304 The fault is defined by a starting point V0 and a list of ``strikes`` and length ``ls``. The strikes and the length
305 is used to define a polyline with points ``V[i]`` such that
306
307 - ``V[0]=V0``
308 - ``V[i]=V[i]+ls[i]*array(cos(strikes[i]),sin(strikes[i]),0)''
309
310 So ``strikes`` defines the angle between the direction of the fault segment and the x0 axis. In 3D ''ls[i]``==0 is allowed.
311
312 In case of a 3D model a fault plane is defined through a dip ``dips`` and depth ``depths``.
313 From the dip and the depth the polyline ``bottom`` of the bottom of the fault is computed.
314
315
316 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]``
317 The segment is parametrized by ``w0`` and ``w1`` with ``w0_offsets[i]<=w0<=w0_offsets[i+1]`` and ``-w1_max<=w1<=0``. Moreover
318
319 - ``(w0,w1)=(w0_offsets[i] , 0)->v0``
320 - ``(w0,w1)=(w0_offsets[i+1], 0)->v1``
321 - ``(w0,w1)=(w0_offsets[i] , -w1_max)->v2``
322 - ``(w0,w1)=(w0_offsets[i+1], -w1_max)->v3``
323
324 If no ``w0_offsets`` is given,
325
326 - ``w0_offsets[0]=0``
327 - ``w0_offsets[i]=w0_offsets[i-1]+L[i]``
328
329 where ``L[i]`` is the length of the segments on the top in 2D and in the middle of the segment in 3D.
330
331 If no ``w1_max`` is given, the average fault depth is used.
332
333
334 :param strikes: list of strikes. This is the angle of the fault segment direction with x0 axis. Right hand rule applies.
335 :type strikes: ``list`` of ``float``
336 :param ls: list of fault lengths. In the case of a 3D fault a segment may have length 0.
337 :type ls: ``list`` of ``float``
338 :param V0: start point of the fault
339 :type V0: ``list`` or `numpy.ndarray' with 2 or 3 components. ``V0[2]`` must be zero.
340 :param tag: the tag of the fault. If fault ``tag`` already exists it is overwritten.
341 :type tag: ``float`` or ``str``
342 :param dips: list of dip angles. Right hand rule around strike direction applies.
343 :type dips: ``list`` of ``float``
344 :param depths: list of segment depth. Value mut be positive in the 3D case.
345 :type depths: ``list`` of ``float``
346 :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.
347 :type w0_offsets: ``list`` of ``float`` or ``None``
348 :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.
349 :type w1_max: ``float``
350 :note: In the three dimensional case the lists ``dip`` and ``top`` must have the same length.
351 """
352 if tag==None:
353 tag=self.NOTAG
354 else:
355 if self.NOTAG in self.getTags():
356 raise ValueError,'Attempt to add a fault with no tag to a set of existing faults'
357 if not isinstance(strikes, list): strikes=[strikes, ]
358 n_segs=len(strikes)
359 if not isinstance(ls, list): ls=[ ls for i in xrange(n_segs) ]
360 if not n_segs==len(ls):
361 raise ValueError,"number of strike direction and length must match."
362 if len(V0)>2:
363 if abs(V0[2])>0: raise Value,"start point needs to be surface (3rd component ==0)"
364 if self.getDim()==2 and not (dips==None and depths == None) :
365 raise ValueError,'Spatial dimension two does not support dip and depth for faults.'
366 if not dips == None:
367 if not isinstance(dips, list): dips=[dips for i in xrange(n_segs) ]
368 if n_segs != len(dips):
369 raise ValueError,'length of dips must be one less then the length of top.'
370 if not depths == None:
371 if not isinstance(depths, list): depths=[depths for i in xrange(n_segs+1) ]
372 if n_segs+1 != len(depths):
373 raise ValueError,'length of depths must be one less then the length of top.'
374 if w0_offsets != None:
375 if len(w0_offsets) != n_segs+1:
376 raise ValueError,'expected length of w0_offsets is %s'%(n_segs)
377 self.__center=None
378 self.__orientation = None
379 #
380 # in the 2D case we don't allow zero length:
381 #
382 if self.getDim() == 2:
383 for l in ls:
384 if l<=0: raise ValueError,"length must be positive"
385 else:
386 for l in ls:
387 if l<0: raise ValueError,"length must be non-negative"
388 for i in xrange(n_segs+1):
389 if depths[i]<0: raise ValueError,"negative depth."
390 #
391 # translate start point to numarray
392 #
393 V0= numpy.array(V0[:self.getDim()],numpy.double)
394 #
395 # set strike vectors:
396 #
397 strike_vectors=[]
398 top_polyline=[V0]
399 total_length=0
400 for i in xrange(n_segs):
401 v=numpy.zeros((self.getDim(),))
402 v[0]=cos(strikes[i])
403 v[1]=sin(strikes[i])
404 strike_vectors.append(v)
405 top_polyline.append(top_polyline[-1]+ls[i]*v)
406 total_length+=ls[i]
407 #
408 # normal and depth direction
409 #
410 if self.getDim()==3:
411 normals=[]
412 for i in xrange(n_segs):
413 normals.append(numpy.array([sin(dips[i])*strike_vectors[i][1], -sin(dips[i])*strike_vectors[i][0], cos(dips[i])]) )
414
415 d=numpy.cross(strike_vectors[0],normals[0])
416 if d[2]>0:
417 f=-1
418 else:
419 f=1
420 depth_vectors=[f*depths[0]*d/numpy.linalg.norm(d) ]
421 for i in xrange(1,n_segs):
422 d=-numpy.cross(normals[i-1],normals[i])
423 d_l=numpy.linalg.norm(d)
424 if d_l<=0:
425 d=numpy.cross(strike_vectors[i],normals[i])
426 d_l=numpy.linalg.norm(d)
427 else:
428 for L in [ strike_vectors[i], strike_vectors[i-1]]:
429 if numpy.linalg.norm(numpy.cross(L,d)) <= self.MIN_DEPTH_ANGLE * numpy.linalg.norm(L) * d_l:
430 raise ValueError,"%s-th depth vector %s too flat."%(i, d)
431 if d[2]>0:
432 f=-1
433 else:
434 f=1
435 depth_vectors.append(f*d*depths[i]/d_l)
436 d=numpy.cross(strike_vectors[n_segs-1],normals[n_segs-1])
437 if d[2]>0:
438 f=-1
439 else:
440 f=1
441 depth_vectors.append(f*depths[n_segs]*d/numpy.linalg.norm(d))
442 bottom_polyline=[ top_polyline[i]+depth_vectors[i] for i in xrange(n_segs+1) ]
443 #
444 # calculate offsets if required:
445 #
446 if w0_offsets==None:
447 w0_offsets=[0.]
448 for i in xrange(n_segs):
449 if self.getDim()==3:
450 w0_offsets.append(w0_offsets[-1]+(float(numpy.linalg.norm(bottom_polyline[i+1]-bottom_polyline[i]))+ls[i])/2.)
451 else:
452 w0_offsets.append(w0_offsets[-1]+ls[i])
453 w0_max=max(w0_offsets)
454 if self.getDim()==3:
455 self.__normals[tag]=normals
456 self.__depth_vectors[tag]=depth_vectors
457 self.__depths[tag]=depths
458 self.__dips[tag]=dips
459 self.__bottom[tag]=bottom_polyline
460 self.__ls[tag]=ls
461 self.__strikes[tag]=strikes
462 self.__strike_vectors[tag]=strike_vectors
463 self.__top[tag]=top_polyline
464 self.__total_length[tag]=total_length
465 self.__offsets[tag]=w0_offsets
466
467 if self.getDim()==2:
468 self.__medDepth[tag]=0.
469 else:
470 self.__medDepth[tag]=sum([ numpy.linalg.norm(v) for v in depth_vectors])/len(depth_vectors)
471 if w1_max==None or self.getDim()==2: w1_max=self.__medDepth[tag]
472 self.__w0_max[tag]=w0_max
473 self.__w1_max[tag]=w1_max
474
475 def getMaxValue(self,f, tol=sqrt(EPSILON)):
476 """
477 returns the maximum value of ``f``, the fault and the location on the fault in fault coordinates.
478
479 :param f: a distribution of values
480 :type f: `escript.Data`
481 :param tol: relative tolerance used to decide if point is on fault
482 :type f: ``tol``
483 :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
484 """
485 ref=-Lsup(f)*2
486 f_max=ref
487 t_max=None
488 p_max=None
489 x=f.getFunctionSpace().getX()
490 for t in self.getTags():
491 p,m=self.getParametrization(x,tag=t, tol=tol)
492 loc=((m*f)+(1.-m)*ref).maxGlobalDataPoint()
493 f_t=f.getTupleForGlobalDataPoint(*loc)[0]
494 if f_t>f_max:
495 f_max=f_t
496 t_max=t
497 p_max=p.getTupleForGlobalDataPoint(*loc)[0]
498
499 return f_max, t_max, p_max
500
501 def getMinValue(self,f, tol=sqrt(EPSILON)):
502 """
503 returns the minimum value of ``f``, the fault and the location on the fault in fault coordinates.
504
505 :param f: a distribution of values
506 :type f: `escript.Data`
507 :param tol: relative tolerance used to decide if point is on fault
508 :type f: ``tol``
509 :return: the minimum value of across all faults in the fault system, the fault tag the minimum 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
510 """
511 ref=Lsup(f)*2
512 f_min=ref
513 t_min=None
514 p_min=None
515 x=f.getFunctionSpace().getX()
516 for t in self.getTags():
517 p,m=self.getParametrization(x,tag=t, tol=tol)
518 loc=((m*f)+(1.-m)*ref).minGlobalDataPoint()
519 f_t=f.getTupleForGlobalDataPoint(*loc)[0]
520 if f_t<f_min:
521 f_min=f_t
522 t_min=t
523 p_min=p.getTupleForGlobalDataPoint(*loc)[0]
524
525 return f_min, t_min, p_min
526
527 def getParametrization(self,x,tag=None, tol=sqrt(EPSILON), outsider=None):
528 """
529 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``.
530
531 Typical usage of the this method is
532
533 dom=Domain(..)
534 x=dom.getX()
535 fs=FaultSystem()
536 fs.addFault(tag=3,...)
537 p, m=fs.getParametrization(x, outsider=0,tag=3)
538 saveDataCSV('x.csv',p=p, x=x, mask=m)
539
540 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.
541
542 :param x: location(s)
543 :type x: `escript.Data` object or `numpy.ndarray`
544 :param tag: the tag of the fault
545 :param tol: relative tolerance to check if location is on fault.
546 :type tol: ``float``
547 :param outsider: value used for parametrization values outside the fault. If not present an appropriate value is choosen.
548 :type outsider: ``float``
549 :return: the coordinates ``x`` in the coordinate system of the fault and a mask indicating coordinates in the fault by 1 (0 elsewhere)
550 :rtype: `escript.Data` object or `numpy.ndarray`
551 """
552 offsets=self.getW0Offsets(tag)
553 w1_range=self.getW1Range(tag)
554 w0_range=self.getW0Range(tag)[1]-self.getW0Range(tag)[0]
555 if outsider == None:
556 outsider=min(self.getW0Range(tag)[0],self.getW0Range(tag)[1])-abs(w0_range)/sqrt(EPSILON)
557
558 if isinstance(x,list): x=numpy.array(x, numpy.double)
559 updated=x[0]*0
560
561 if self.getDim()==2:
562 #
563 #
564 p=x[0]*0 + outsider
565 top=self.getTopPolyline(tag)
566 for i in xrange(1,len(top)):
567 d=top[i]-top[i-1]
568 h=x-top[i-1]
569 h_l=length(h)
570 d_l=length(d)
571 s=inner(h,d)/d_l**2
572 s=s*whereNonPositive(s-1.-tol)*whereNonNegative(s+tol)
573 m=whereNonPositive(length(h-s*d)-tol*maximum(h_l,d_l))*(1.-updated)
574 p=(1.-m)*p+m*(offsets[i-1]+(offsets[i]-offsets[i-1])*s)
575 updated=wherePositive(updated+m)
576 else:
577 p=x[:2]*0 + outsider
578 top=self.getTopPolyline(tag)
579 bottom=self.getBottomPolyline(tag)
580 n=self.getSegmentNormals(tag)
581 for i in xrange(len(top)-1):
582 h=x-top[i]
583 R=top[i+1]-top[i]
584 r=bottom[i+1]-bottom[i]
585 D0=bottom[i]-top[i]
586 D1=bottom[i+1]-top[i+1]
587 s_upper=matrix_mult(numpy.linalg.pinv(numpy.vstack((R,D1)).T),h)
588 s_lower=matrix_mult(numpy.linalg.pinv(numpy.vstack((r,D0)).T),h)
589 m_ul=wherePositive(s_upper[0]-s_upper[1])
590 s=s_upper*m_ul+s_lower*(1-m_ul)
591 s0=s[0]
592 s1=s[1]
593 m=whereNonNegative(s0+tol)*whereNonPositive(s0-1.-tol)*whereNonNegative(s1+tol)*whereNonPositive(s1-1.-tol)
594 s0=s0*m
595 s1=s1*m
596 atol=tol*maximum(length(h),length(top[i]-bottom[i+1]))
597 m=whereNonPositive(length(h-s0*R-s1*D1)*m_ul+(1-m_ul)*length(h-s0*r-s1*D0)-atol)
598 p[0]=(1.-m)*p[0]+m*(offsets[i]+(offsets[i+1]-offsets[i])*s0)
599 p[1]=(1.-m)*p[1]+m*(w1_range[1]+(w1_range[0]-w1_range[1])*s1)
600 updated=wherePositive(updated+m)
601
602 return p, updated
603
604 def getSideAndDistance(self,x,tag=None):
605 """
606 returns the side and the distance at ``x`` from the fault ``tag``.
607
608 :param x: location(s)
609 :type x: `escript.Data` object or `numpy.ndarray`
610 :param tag: the tag of the fault
611 :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.
612 """
613 d=None
614 side=None
615 if self.getDim()==2:
616 mat=numpy.array([[0., 1.], [-1., 0.] ])
617 s=self.getTopPolyline(tag)
618 for i in xrange(1,len(s)):
619 q=(s[i]-s[i-1])
620 h=x-s[i-1]
621 q_l=length(q)
622 qt=matrixmult(mat,q) # orthogonal direction
623 t=inner(q,h)/q_l**2
624 t=maximum(minimum(t,1,),0.)
625 p=h-t*q
626 dist=length(p)
627 lside=sign(inner(p,qt))
628 if d == None:
629 d=dist
630 side=lside
631 else:
632 m=whereNegative(d-dist)
633 m2=wherePositive(whereZero(abs(lside))+m)
634 d=dist*(1-m)+d*m
635 side=lside*(1-m2)+side*m2
636 else:
637 ns=self.getSegmentNormals(tag)
638 top=self.getTopPolyline(tag)
639 bottom=self.getBottomPolyline(tag)
640 for i in xrange(len(top)-1):
641 h=x-top[i]
642 R=top[i+1]-top[i]
643 r=bottom[i+1]-bottom[i]
644 D0=bottom[i]-top[i]
645 D1=bottom[i+1]-top[i+1]
646 s_upper=matrix_mult(numpy.linalg.pinv(numpy.vstack((R,D1)).T),h)
647 s_lower=matrix_mult(numpy.linalg.pinv(numpy.vstack((r,D0)).T),h)
648 m_ul=wherePositive(s_upper[0]-s_upper[1])
649 s=s_upper*m_ul+s_lower*(1-m_ul)
650 s=maximum(minimum(s,1.),0)
651 p=h-s[0]*(R*m_ul+(1-m_ul)*r)-s[1]*(D1*m_ul+(1-m_ul)*D0)
652 dist=length(p)
653 lside=sign(inner(p,ns[i]))
654 if d == None:
655 d=dist
656 side=lside
657 else:
658 m=whereNegative(d-dist)
659 m2=wherePositive(whereZero(abs(lside))+m)
660 d=dist*(1-m)+d*m
661 side=lside*(1-m2)+side*m2
662
663 return side, d
664

  ViewVC Help
Powered by ViewVC 1.1.26