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

  ViewVC Help
Powered by ViewVC 1.1.26