/[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 2663 - (show annotations)
Tue Sep 15 06:04:44 2009 UTC (10 years, 2 months ago) by gross
File MIME type: text/x-python
File size: 23616 byte(s)
finally there is a 3D version for the fault system class
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 depth_vectors=[depths[0]*d/numpy.linalg.norm(d) ]
417 for i in xrange(1,n_segs):
418 d=-numpy.cross(normals[i-1],normals[i])
419 d_l=numpy.linalg.norm(d)
420 if d_l<=0:
421 d=numpy.cross(strike_vectors[i],normals[i])
422 depth_vectors.append([depths[d]*d/numpy.linalg.norm(d)])
423 else:
424 d=depths[i]*d/d_l
425 for L in [ strike_vectors[i], strike_vectors[i-1]]:
426 if numpy.linalg.norm(numpy.cross(L,d)) <= self.MIN_DEPTH_ANGLE * numpy.linalg.norm(L) * depths[i]:
427 raise ValueError,"%s-th depth vector %s too flat."%(i, d)
428 depth_vectors.append(d)
429 d=numpy.cross(strike_vectors[n_segs-1],normals[n_segs-1])
430 depth_vectors.append(depths[n_segs]*d/numpy.linalg.norm(d))
431 bottom_polyline=[ top_polyline[i]+depth_vectors[i] for i in xrange(n_segs+1) ]
432 #
433 # calculate offsets if required:
434 #
435 if w0_offsets==None:
436 w0_offsets=[0.]
437 for i in xrange(n_segs):
438 if self.getDim()==3:
439 w0_offsets.append(w0_offsets[-1]+(float(numpy.linalg.norm(bottom_polyline[i+1]-bottom_polyline[i]))+ls[i])/2.)
440 else:
441 w0_offsets.append(w0_offsets[-1]+ls[i])
442 w0_max=max(w0_offsets)
443 if self.getDim()==3:
444 self.__normals[tag]=normals
445 self.__depth_vectors[tag]=depth_vectors
446 self.__depths[tag]=depths
447 self.__dips[tag]=dips
448 self.__bottom[tag]=bottom_polyline
449 self.__ls[tag]=ls
450 self.__strikes[tag]=strikes
451 self.__strike_vectors[tag]=strike_vectors
452 self.__top[tag]=top_polyline
453 self.__total_length[tag]=total_length
454 self.__offsets[tag]=w0_offsets
455
456 if self.getDim()==2:
457 self.__medDepth[tag]=0.
458 else:
459 self.__medDepth[tag]=sum([ numpy.linalg.norm(v) for v in depth_vectors])/len(depth_vectors)
460 if w1_max==None or self.getDim()==2: w1_max=self.__medDepth[tag]
461 self.__w0_max[tag]=w0_max
462 self.__w1_max[tag]=w1_max
463
464 def getMaxValue(self,f, tol=sqrt(EPSILON)):
465 """
466 returns the maximum value of ``f``, the fault and the location on the fault in fault coordinates.
467
468 :param f: a distribution of values
469 :type f: `escript.Data`
470 :param tol: relative tolerance used to decide if point is on fault
471 :type f: ``tol``
472 :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
473 """
474 ref=-Lsup(f)*2
475 f_max=ref
476 t_max=None
477 p_max=None
478 x=f.getFunctionSpace().getX()
479 for t in self.getTags():
480 p,m=self.getParametrization(x,tag=t, tol=tol)
481 loc=((m*f)+(1.-m)*ref).maxGlobalDataPoint()
482 f_t=f.getTupleForGlobalDataPoint(*loc)[0]
483 if f_t>f_max:
484 f_max=f_t
485 t_max=t
486 p_max=p.getTupleForGlobalDataPoint(*loc)[0]
487
488 return f_max, t_max, p_max
489
490 def getParametrization(self,x,tag=None, tol=sqrt(EPSILON), outsider=None):
491 """
492 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``.
493
494 Typical usage of the this method is
495
496 dom=Domain(..)
497 x=dom.getX()
498 fs=FaultSystem()
499 fs.addFault(tag=3,...)
500 p, m=fs.getParametrization(x, outsider=0,tag=3)
501 saveDataCSV('x.csv',p=p, x=x, mask=m)
502
503 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.
504
505 :param x: location(s)
506 :type x: `escript.Data` object or `numpy.ndarray`
507 :param tag: the tag of the fault
508 :param tol: relative tolerance to check if location is on fault.
509 :type tol: ``float``
510 :param outsider: value used for parametrization values outside the fault. If not present an appropriate value is choosen.
511 :type outsider: ``float``
512 :return: the coordinates ``x`` in the coordinate system of the fault and a mask indicating coordinates in the fault by 1 (0 elsewhere)
513 :rtype: `escript.Data` object or `numpy.ndarray`
514 """
515 offsets=self.getW0Offsets(tag)
516 w1_range=self.getW1Range(tag)
517 w0_range=self.getW0Range(tag)[1]-self.getW0Range(tag)[0]
518 if outsider == None:
519 outsider=min(self.getW0Range(tag)[0],self.getW0Range(tag)[1])-abs(w0_range)/sqrt(EPSILON)
520
521 if isinstance(x,list): x=numpy.array(x, numpy.double)
522 updated=x[0]*0
523
524 if self.getDim()==2:
525 #
526 #
527 p=x[0]*0 + outsider
528 top=self.getTopPolyline(tag)
529 for i in xrange(1,len(top)):
530 d=top[i]-top[i-1]
531 h=x-top[i-1]
532 h_l=length(h)
533 d_l=length(d)
534 s=inner(h,d)/d_l**2
535 s=s*whereNonPositive(s-1.-tol)*whereNonNegative(s+tol)
536 m=whereNonPositive(length(h-s*d)-tol*maximum(h_l,d_l))*(1.-updated)
537 p=(1.-m)*p+m*(offsets[i-1]+(offsets[i]-offsets[i-1])*s)
538 updated=wherePositive(updated+m)
539 else:
540 p=x[:2]*0 + outsider
541 top=self.getTopPolyline(tag)
542 bottom=self.getBottomPolyline(tag)
543 n=self.getSegmentNormals(tag)
544 for i in xrange(len(top)-1):
545 h=x-top[i]
546 R=top[i+1]-top[i]
547 r=bottom[i+1]-bottom[i]
548 D0=bottom[i]-top[i]
549 D1=bottom[i+1]-top[i+1]
550 s_upper=matrix_mult(numpy.linalg.pinv(numpy.vstack((R,D1)).T),h)
551 s_lower=matrix_mult(numpy.linalg.pinv(numpy.vstack((r,D0)).T),h)
552 m_ul=wherePositive(s_upper[0]-s_upper[1])
553 s=s_upper*m_ul+s_lower*(1-m_ul)
554 s0=s[0]
555 s1=s[1]
556 m=whereNonNegative(s0+tol)*whereNonPositive(s0-1.-tol)*whereNonNegative(s1+tol)*whereNonPositive(s1-1.-tol)
557 s0=s0*m
558 s1=s1*m
559 atol=tol*maximum(length(h),length(top[i]-bottom[i+1]))
560 m=whereNonPositive(length(h-s0*R-s1*D1)*m_ul+(1-m_ul)*length(h-s0*r-s1*D0)-atol)
561 p[0]=(1.-m)*p[0]+m*(offsets[i]+(offsets[i+1]-offsets[i])*s0)
562 p[1]=(1.-m)*p[1]+m*(w1_range[1]+(w1_range[0]-w1_range[1])*s1)
563 updated=wherePositive(updated+m)
564
565 return p, updated
566
567 def getSideAndDistance(self,x,tag=None):
568 """
569 returns the side and the distance at ``x`` from the fault ``tag``.
570
571 :param x: location(s)
572 :type x: `escript.Data` object or `numpy.ndarray`
573 :param tag: the tag of the fault
574 :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.
575 """
576 d=None
577 side=None
578 if self.getDim()==2:
579 mat=numpy.array([[0., 1.], [-1., 0.] ])
580 s=self.getTopPolyline(tag)
581 for i in xrange(1,len(s)):
582 q=(s[i]-s[i-1])
583 h=x-s[i-1]
584 q_l=length(q)
585 qt=matrixmult(mat,q) # orthogonal direction
586 t=inner(q,h)/q_l**2
587 t=maximum(minimum(t,1,),0.)
588 p=h-t*q
589 dist=length(p)
590 lside=sign(inner(p,qt))
591 if d == None:
592 d=dist
593 side=lside
594 else:
595 m=whereNegative(d-dist)
596 m2=wherePositive(whereZero(abs(lside))+m)
597 d=dist*(1-m)+d*m
598 side=lside*(1-m2)+side*m2
599 else:
600 ns=self.getSegmentNormals(tag)
601 top=self.getTopPolyline(tag)
602 bottom=self.getBottomPolyline(tag)
603 for i in xrange(len(top)-1):
604 h=x-top[i]
605 R=top[i+1]-top[i]
606 r=bottom[i+1]-bottom[i]
607 D0=bottom[i]-top[i]
608 D1=bottom[i+1]-top[i+1]
609 s_upper=matrix_mult(numpy.linalg.pinv(numpy.vstack((R,D1)).T),h)
610 s_lower=matrix_mult(numpy.linalg.pinv(numpy.vstack((r,D0)).T),h)
611 m_ul=wherePositive(s_upper[0]-s_upper[1])
612 s=s_upper*m_ul+s_lower*(1-m_ul)
613 s=maximum(minimum(s,1.),0)
614 p=h-s[0]*(R*m_ul+(1-m_ul)*r)-s[1]*(D1*m_ul+(1-m_ul)*D0)
615 dist=length(p)
616 lside=sign(inner(p,ns[i]))
617 if d == None:
618 d=dist
619 side=lside
620 else:
621 m=whereNegative(d-dist)
622 m2=wherePositive(whereZero(abs(lside))+m)
623 d=dist*(1-m)+d*m
624 side=lside*(1-m2)+side*m2
625
626 return side, d
627

  ViewVC Help
Powered by ViewVC 1.1.26