/[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 2675 - (hide 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 gross 2647 ########################################################
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 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
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 gross 2647 """
41     NOTAG="__NOTAG__"
42 gross 2663 MIN_DEPTH_ANGLE=0.1
43 gross 2647 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 gross 2663 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 gross 2647 self.__offsets={}
72     self.__w1_max={}
73     self.__w0_max={}
74     self.__center=None
75     self.__orientation = None
76 gross 2663 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 gross 2647
83 gross 2663 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 gross 2647 def getDim(self):
90     """
91     returns the spatial dimension
92     :rtype: ``int``
93     """
94     return self.__dim
95 gross 2663
96     def getTopPolyline(self, tag=None):
97 gross 2647 """
98 gross 2663 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 gross 2647 """
104     if tag==None: tag=self.NOTAG
105 gross 2663 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 gross 2647
128 gross 2663 def getTotalLength(self, tag=None):
129 gross 2647 """
130 gross 2663 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 gross 2647 returns the medium depth of fault ``tag``
139     :rtype: ``float``
140     """
141     if tag==None: tag=self.NOTAG
142 gross 2663 return self.__medDepth[tag]
143 gross 2647
144 gross 2663 def getDips(self, tag=None):
145 gross 2647 """
146 gross 2663 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 gross 2647 """
151 gross 2663 if tag==None: tag=self.NOTAG
152     if self.getDim()==3:
153     return self.__dips[tag]
154     else:
155     return None
156 gross 2647
157 gross 2663 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 gross 2647 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 gross 2663 for s in self.getTopPolyline(t):
244 gross 2647 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 gross 2663 for s in self.getTopPolyline(t):
259 gross 2647 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 gross 2663 d=V[:,0]
266     else:
267 gross 2647 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 gross 2663
287 gross 2647 for t in self.getTags():
288 gross 2663 strikes=[ s+ rot for s in self.getStrikes(t) ]
289     V0=self.getStart(t)
290 gross 2647
291 gross 2663 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 gross 2647 """
302 gross 2663 adds a new fault to the fault system. The fault is named by ``tag``.
303 gross 2647
304 gross 2663 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 gross 2647
319     - ``(w0,w1)=(w0_offsets[i] , 0)->v0``
320 gross 2663 - ``(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 gross 2647
324     If no ``w0_offsets`` is given,
325    
326     - ``w0_offsets[0]=0``
327 gross 2663 - ``w0_offsets[i]=w0_offsets[i-1]+L[i]``
328 gross 2647
329 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.
330 gross 2647
331 gross 2663 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 gross 2647 :param tag: the tag of the fault. If fault ``tag`` already exists it is overwritten.
341     :type tag: ``float`` or ``str``
342 gross 2663 :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 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.
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 gross 2663 :note: In the three dimensional case the lists ``dip`` and ``top`` must have the same length.
351 gross 2647 """
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 gross 2663 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 gross 2647 if w0_offsets != None:
375 gross 2663 if len(w0_offsets) != n_segs+1:
376     raise ValueError,'expected length of w0_offsets is %s'%(n_segs)
377 gross 2647 self.__center=None
378     self.__orientation = None
379     #
380 gross 2663 # 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 gross 2647 else:
386 gross 2663 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 gross 2647 #
391 gross 2663 # translate start point to numarray
392 gross 2647 #
393 gross 2663 V0= numpy.array(V0[:self.getDim()],numpy.double)
394 gross 2647 #
395 gross 2663 # set strike vectors:
396 gross 2647 #
397 gross 2663 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 gross 2665 if d[2]>0:
417     f=-1
418     else:
419     f=1
420     depth_vectors=[f*depths[0]*d/numpy.linalg.norm(d) ]
421 gross 2663 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 gross 2665 d_l=numpy.linalg.norm(d)
427 gross 2663 else:
428     for L in [ strike_vectors[i], strike_vectors[i-1]]:
429 gross 2665 if numpy.linalg.norm(numpy.cross(L,d)) <= self.MIN_DEPTH_ANGLE * numpy.linalg.norm(L) * d_l:
430 gross 2663 raise ValueError,"%s-th depth vector %s too flat."%(i, d)
431 gross 2665 if d[2]>0:
432     f=-1
433     else:
434     f=1
435     depth_vectors.append(f*d*depths[i]/d_l)
436 gross 2663 d=numpy.cross(strike_vectors[n_segs-1],normals[n_segs-1])
437 gross 2665 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 gross 2663 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 gross 2647 else:
470 gross 2663 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 gross 2647 self.__w0_max[tag]=w0_max
473 gross 2663 self.__w1_max[tag]=w1_max
474 gross 2647
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 gross 2675 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 gross 2647 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 gross 2663 :param tag: the tag of the fault
545 gross 2647 :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 gross 2663 w1_range=self.getW1Range(tag)
554 gross 2647 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 gross 2663 if isinstance(x,list): x=numpy.array(x, numpy.double)
559 gross 2647 updated=x[0]*0
560    
561     if self.getDim()==2:
562     #
563     #
564 gross 2663 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 gross 2647 h_l=length(h)
570     d_l=length(d)
571 gross 2663 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 gross 2647 updated=wherePositive(updated+m)
576     else:
577 gross 2663 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 gross 2647
602     return p, updated
603 gross 2649
604 gross 2654 def getSideAndDistance(self,x,tag=None):
605 gross 2649 """
606 gross 2654 returns the side and the distance at ``x`` from the fault ``tag``.
607 gross 2647
608 gross 2649 :param x: location(s)
609     :type x: `escript.Data` object or `numpy.ndarray`
610 gross 2663 :param tag: the tag of the fault
611 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.
612 gross 2649 """
613 gross 2654 d=None
614     side=None
615 gross 2649 if self.getDim()==2:
616 gross 2654 mat=numpy.array([[0., 1.], [-1., 0.] ])
617 gross 2663 s=self.getTopPolyline(tag)
618 gross 2649 for i in xrange(1,len(s)):
619 gross 2654 q=(s[i]-s[i-1])
620 gross 2649 h=x-s[i-1]
621 gross 2654 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 gross 2649 else:
637 gross 2663 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 gross 2654 return side, d
664 gross 2647

  ViewVC Help
Powered by ViewVC 1.1.26