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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 3360 - (hide annotations)
Thu Nov 18 00:20:21 2010 UTC (9 years, 11 months ago) by jfenwick
File MIME type: text/x-python
File size: 24968 byte(s)
Fix some epydoc errors
Fixed issue 564
removed minimize.py since it hasn't worked in who knows how long
1 gross 2647 ########################################################
2     #
3 jfenwick 2881 # Copyright (c) 2003-2010 by University of Queensland
4 gross 2647 # 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 jfenwick 2881 __copyright__="""Copyright (c) 2003-2010 by University of Queensland
14 gross 2647 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 gross 2676 from esys.escript.pdetools import Locator
23 gross 2647 import numpy
24     import math
25    
26     class FaultSystem:
27     """
28 jfenwick 3360 The FaultSystem class defines a system of faults in the Earth's crust.
29    
30 gross 2663 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 jfenwick 3360 - ``V[i]=V[i]+ls[i]*array(cos(strike[i]),sin(strike[i]),0)``
35 gross 2663
36 jfenwick 3360 So ``strike`` defines the angle between the direction of the fault segment and the x0 axis. ls[i]==0 is allowed.
37 gross 2663
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 gross 2647 """
42     NOTAG="__NOTAG__"
43 gross 2663 MIN_DEPTH_ANGLE=0.1
44 gross 2647 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 gross 2663 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 gross 2647 self.__offsets={}
73     self.__w1_max={}
74     self.__w0_max={}
75     self.__center=None
76     self.__orientation = None
77 gross 2663 def getStart(self,tag=None):
78     """
79     returns the starting point of fault ``tag``
80 jfenwick 3360 :rtype: `numpy.ndarray`.
81 gross 2663 """
82     return self.getTopPolyline(tag)[0]
83 gross 2647
84 gross 2663 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 gross 2647 def getDim(self):
91     """
92     returns the spatial dimension
93     :rtype: ``int``
94     """
95     return self.__dim
96 gross 2663
97     def getTopPolyline(self, tag=None):
98 gross 2647 """
99 gross 2663 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 jfenwick 3360 :return: the list of vertices defining the top of the fault. The coordinates are `numpy.ndarray`.
104 gross 2647 """
105     if tag==None: tag=self.NOTAG
106 gross 2663 return self.__top[tag]
107     def getStrikes(self, tag=None):
108     """
109 jfenwick 3360 :return: the strike of the segements in fault ``tag``
110 gross 2663 :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 jfenwick 3360 :return: the strike vectors of fault ``tag``
117     :rtype: ``list`` of `numpy.ndarray`.
118 gross 2663 """
119     if tag==None: tag=self.NOTAG
120     return self.__strike_vectors[tag]
121     def getLengths(self, tag=None):
122     """
123 jfenwick 3360 :return: the lengths of segments in fault ``tag``
124     :rtype: ``list`` of ``float``
125 gross 2663 """
126     if tag==None: tag=self.NOTAG
127     return self.__ls[tag]
128 gross 2647
129 gross 2663 def getTotalLength(self, tag=None):
130 gross 2647 """
131 jfenwick 3360 :return: the total unrolled length of fault ``tag``
132     :rtype: ``float``
133 gross 2663 """
134     if tag==None: tag=self.NOTAG
135     return self.__total_length[tag]
136    
137     def getMediumDepth(self,tag=None):
138     """
139 gross 2647 returns the medium depth of fault ``tag``
140     :rtype: ``float``
141     """
142     if tag==None: tag=self.NOTAG
143 gross 2663 return self.__medDepth[tag]
144 gross 2647
145 gross 2663 def getDips(self, tag=None):
146 gross 2647 """
147 gross 2663 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 gross 2647 """
152 gross 2663 if tag==None: tag=self.NOTAG
153     if self.getDim()==3:
154     return self.__dips[tag]
155     else:
156     return None
157 gross 2647
158 gross 2663 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 gross 2647 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 gross 2663 for s in self.getTopPolyline(t):
245 gross 2647 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 gross 2663 for s in self.getTopPolyline(t):
260 gross 2647 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 gross 2663 d=V[:,0]
267     else:
268 gross 2647 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 gross 2663
288 gross 2647 for t in self.getTags():
289 gross 2663 strikes=[ s+ rot for s in self.getStrikes(t) ]
290     V0=self.getStart(t)
291 gross 2647
292 gross 2663 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 gross 2647 """
303 gross 2663 adds a new fault to the fault system. The fault is named by ``tag``.
304 gross 2647
305 gross 2663 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 jfenwick 3360 - ``V[i]=V[i]+ls[i]*array(cos(strikes[i]),sin(strikes[i]),0)``
310 gross 2663
311 jfenwick 3360 So ``strikes`` defines the angle between the direction of the fault segment and the x0 axis. In 3D ``ls[i]``==0 is allowed.
312 gross 2663
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 jfenwick 3360 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 gross 2663 The segment is parametrized by ``w0`` and ``w1`` with ``w0_offsets[i]<=w0<=w0_offsets[i+1]`` and ``-w1_max<=w1<=0``. Moreover
319 gross 2647
320     - ``(w0,w1)=(w0_offsets[i] , 0)->v0``
321 gross 2663 - ``(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 gross 2647
325     If no ``w0_offsets`` is given,
326    
327     - ``w0_offsets[0]=0``
328 gross 2663 - ``w0_offsets[i]=w0_offsets[i-1]+L[i]``
329 gross 2647
330 gross 2663 where ``L[i]`` is the length of the segments on the top in 2D and in the middle of the segment in 3D.
331 gross 2647
332 gross 2663 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 gross 2647 :param tag: the tag of the fault. If fault ``tag`` already exists it is overwritten.
342     :type tag: ``float`` or ``str``
343 gross 2663 :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 gross 2647 :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 gross 2663 :note: In the three dimensional case the lists ``dip`` and ``top`` must have the same length.
352 gross 2647 """
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 gross 2663 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 gross 2647 if w0_offsets != None:
376 gross 2663 if len(w0_offsets) != n_segs+1:
377     raise ValueError,'expected length of w0_offsets is %s'%(n_segs)
378 gross 2647 self.__center=None
379     self.__orientation = None
380     #
381 gross 2663 # 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 gross 2647 else:
387 gross 2663 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 gross 2647 #
392 gross 2663 # translate start point to numarray
393 gross 2647 #
394 gross 2663 V0= numpy.array(V0[:self.getDim()],numpy.double)
395 gross 2647 #
396 gross 2663 # set strike vectors:
397 gross 2647 #
398 gross 2663 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 gross 2676 normals.append(numpy.array([sin(dips[i])*strike_vectors[i][1],-sin(dips[i])*strike_vectors[i][0], cos(dips[i])]) )
415 gross 2663
416     d=numpy.cross(strike_vectors[0],normals[0])
417 gross 2665 if d[2]>0:
418     f=-1
419     else:
420     f=1
421     depth_vectors=[f*depths[0]*d/numpy.linalg.norm(d) ]
422 gross 2663 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 gross 2665 d_l=numpy.linalg.norm(d)
428 gross 2663 else:
429     for L in [ strike_vectors[i], strike_vectors[i-1]]:
430 gross 2665 if numpy.linalg.norm(numpy.cross(L,d)) <= self.MIN_DEPTH_ANGLE * numpy.linalg.norm(L) * d_l:
431 gross 2663 raise ValueError,"%s-th depth vector %s too flat."%(i, d)
432 gross 2665 if d[2]>0:
433     f=-1
434     else:
435     f=1
436     depth_vectors.append(f*d*depths[i]/d_l)
437 gross 2663 d=numpy.cross(strike_vectors[n_segs-1],normals[n_segs-1])
438 gross 2665 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 gross 2663 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 gross 2647 else:
471 gross 2663 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 gross 2647 self.__w0_max[tag]=w0_max
474 gross 2663 self.__w1_max[tag]=w1_max
475 gross 2647
476     def getMaxValue(self,f, tol=sqrt(EPSILON)):
477     """
478 gross 2676 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 gross 2647
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 jfenwick 3360 :type tol: ``tol``
484 gross 2676 :return: the fault tag the maximum is taken, and a `Locator` object to collect the value at location of maximum value.
485 gross 2647 """
486     ref=-Lsup(f)*2
487     f_max=ref
488     t_max=None
489 gross 2676 loc_max=None
490 gross 2647 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 gross 2676 loc_max=loc
499 gross 2647
500 gross 2676 if loc_max == None:
501     return None, None
502     else:
503     return t_max, Locator(x.getFunctionSpace(),x.getTupleForGlobalDataPoint(*loc_max))
504 gross 2647
505 gross 2675 def getMinValue(self,f, tol=sqrt(EPSILON)):
506     """
507 gross 2676 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 gross 2675
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 jfenwick 3360 :type tol: ``tol``
513 gross 2676 :return: the fault tag the minimum is taken, and a `Locator` object to collect the value at location of minimum value.
514 gross 2675 """
515     ref=Lsup(f)*2
516     f_min=ref
517     t_min=None
518 gross 2676 loc_min=None
519 gross 2675 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 gross 2676 loc_min=loc
528 gross 2675
529 gross 2676 if loc_min == None:
530     return None, None
531     else:
532     return t_min, Locator(x.getFunctionSpace(),x.getTupleForGlobalDataPoint(*loc_min))
533 gross 2675
534 gross 2647 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 gross 2663 :param tag: the tag of the fault
552 gross 2647 :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 gross 2663 w1_range=self.getW1Range(tag)
561 gross 2647 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 gross 2663 if isinstance(x,list): x=numpy.array(x, numpy.double)
566 gross 2647 updated=x[0]*0
567    
568     if self.getDim()==2:
569     #
570     #
571 gross 2663 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 gross 2647 h_l=length(h)
577     d_l=length(d)
578 gross 2663 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 gross 2647 updated=wherePositive(updated+m)
583     else:
584 gross 2663 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 gross 2647
609     return p, updated
610 gross 2649
611 gross 2654 def getSideAndDistance(self,x,tag=None):
612 gross 2649 """
613 gross 2654 returns the side and the distance at ``x`` from the fault ``tag``.
614 gross 2647
615 gross 2649 :param x: location(s)
616     :type x: `escript.Data` object or `numpy.ndarray`
617 gross 2663 :param tag: the tag of the fault
618 gross 2654 :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 gross 2649 """
620 gross 2654 d=None
621     side=None
622 gross 2649 if self.getDim()==2:
623 gross 2654 mat=numpy.array([[0., 1.], [-1., 0.] ])
624 gross 2663 s=self.getTopPolyline(tag)
625 gross 2649 for i in xrange(1,len(s)):
626 gross 2654 q=(s[i]-s[i-1])
627 gross 2649 h=x-s[i-1]
628 gross 2654 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 gross 2649 else:
644 gross 2663 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 gross 2683 p=h-(m_ul*R+(1-m_ul)*r)*s[0]-(m_ul*D1+(1-m_ul)*D0)*s[1]
659 gross 2663 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 gross 2654 return side, d
671 gross 2647

  ViewVC Help
Powered by ViewVC 1.1.26