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

Contents of /trunk/escriptcore/py_src/faultsystems.py

Parent Directory Parent Directory | Revision Log Revision Log


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

  ViewVC Help
Powered by ViewVC 1.1.26