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

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

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

revision 2662 by gross, Tue Sep 8 07:11:12 2009 UTC revision 2663 by gross, Tue Sep 15 06:04:44 2009 UTC
# Line 26  class FaultSystem: Line 26  class FaultSystem:
26    """    """
27    The FaultSystem class defines a system of faults in the Earth crust.    The FaultSystem class defines a system of faults in the Earth crust.
28    
29    A fault system is defined by set of faults index by a tag. Each fault is defined as a polygon describing the top of    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    the fault and - in case of a 3D model - a polygon describing the bottom of the fault. This class provides a mechanism    strikes ``strike`` and length ``l``. The strikes and the length is used to define a polyline with points ``V[i]`` such that
31    to parametrise each fault with the domain [0,w0_max] x [0, w1_max]  (to [0,w0_max] in the 2D case).  
32      - ``V[0]=V0``
33      - ``V[i]=V[i]+ls[i]*array(cos(strike[i]),sin(strike[i]),0)''
34    
35      So ``strike`` defines the angle between the direction of the fault segment and the x0 axis. ''ls[i]``==0 is allowed.
36    
37      In case of a 3D model a fault plane is defined through a dip and depth.
38    
39      The class provides a mechanism to parametrise each fault with the domain [0,w0_max] x [0, w1_max]  (to [0,w0_max] in the 2D case).
40    """    """
41    NOTAG="__NOTAG__"    NOTAG="__NOTAG__"
42      MIN_DEPTH_ANGLE=0.1
43    def __init__(self,dim=3):    def __init__(self,dim=3):
44      """      """
45      Sets up the fault system      Sets up the fault system
# Line 41  class FaultSystem: Line 50  class FaultSystem:
50      if not (dim == 2 or dim == 3):      if not (dim == 2 or dim == 3):
51         raise ValueError,"only dimension2 2 and 3 are supported."         raise ValueError,"only dimension2 2 and 3 are supported."
52      self.__dim=dim      self.__dim=dim
53      self.__faults={}      self.__top={}
54      self.__length={}      self.__ls={}
55      self.__depth={}      self.__strikes={}
56        self.__strike_vectors={}
57        self.__medDepth={}
58        self.__total_length={}
59        if dim ==2:
60           self.__depths=None
61           self.__depth_vectors=None
62           self.__dips=None
63           self.__bottom=None
64           self.__normals=None
65        else:
66           self.__depths={}
67           self.__depth_vectors={}
68           self.__dips={}
69           self.__bottom={}
70           self.__normals={}
71      self.__offsets={}      self.__offsets={}
72      self.__w1_max={}      self.__w1_max={}
73      self.__w0_max={}      self.__w0_max={}
74      self.__center=None      self.__center=None
75      self.__orientation = None      self.__orientation = None
76      def getStart(self,tag=None):
77         """
78         returns the starting point of fault ``tag``
79         :rtype: `numpy.ndarray'.
80         """
81         return self.getTopPolyline(tag)[0]
82    
83      def getTags(self):
84         """
85         returns a list of the tags used by the fault system
86         :rtype: ``list``
87         """
88         return self.__top.keys()
89    def getDim(self):    def getDim(self):
90       """       """
91       returns the spatial dimension       returns the spatial dimension
92       :rtype: ``int``       :rtype: ``int``
93       """       """
94       return self.__dim       return self.__dim
95    def getLength(self,tag=None):  
96      def getTopPolyline(self, tag=None):
97       """       """
98       returns the unrolled length of fault ``tag``       returns the polyline used to describe fault tagged by ``tag``
99       :rtype: ``float``        
100         :param tag: the tag of the fault
101         :type tag: ``float`` or ``str``
102         :return: the list of vertices defining the top of the fault.  The coordinates are `numpy.ndarray'.
103         """
104         if tag==None: tag=self.NOTAG
105         return self.__top[tag]
106      def getStrikes(self, tag=None):
107         """
108         returns the strike of the segements in fault ``tag``
109         :rtype: ``list`` of ``float``
110         """
111         if tag==None: tag=self.NOTAG
112         return self.__strikes[tag]
113      def getStrikeVectors(self, tag=None):
114         """
115         returns the strike vectors of fault ``tag``
116         :rtype: ``list`` of `numpy.ndarray'.
117         """
118         if tag==None: tag=self.NOTAG
119         return self.__strike_vectors[tag]
120      def getLengths(self, tag=None):
121         """
122         returns the lengths of segments in fault ``tag``
123         :rtype: ``list`` of `float``
124         """
125         if tag==None: tag=self.NOTAG
126         return self.__ls[tag]
127    
128      def getTotalLength(self, tag=None):
129         """
130         returns the total unrolled length of fault ``tag``
131         :rtype: `float``
132       """       """
133       if tag==None: tag=self.NOTAG       if tag==None: tag=self.NOTAG
134       return self.__length[tag]       return self.__total_length[tag]
135    
136    def getDepth(self,tag=None):    def getMediumDepth(self,tag=None):
137       """       """
138       returns the medium depth of fault ``tag``       returns the medium depth of fault ``tag``
139       :rtype: ``float``       :rtype: ``float``
140       """       """
141       if tag==None: tag=self.NOTAG       if tag==None: tag=self.NOTAG
142       return self.__depth[tag]       return self.__medDepth[tag]
143    
144    def getTags(self):    def getDips(self, tag=None):
145       """       """
146       returns a list of the tags used by the fault system       returns the list of the dips of the segements in fault ``tag``
147       :rtype: ``list``       :param tag: the tag of the fault
148         :type tag: ``float`` or ``str``
149         :return: the list of segment dips. In the 2D case None is returned.
150       """       """
151       return self.__faults.keys()       if tag==None: tag=self.NOTAG
152         if self.getDim()==3:
153             return self.__dips[tag]
154         else:
155             return None
156    
157      def getBottomPolyline(self, tag=None):
158         """
159         returns the list of the vertices defining the bottom of the fault ``tag``
160         :param tag: the tag of the fault
161         :type tag: ``float`` or ``str``
162         :return: the list of vertices. In the 2D case None is returned.
163         """
164         if tag==None: tag=self.NOTAG
165         if self.getDim()==3:
166             return self.__bottom[tag]
167         else:
168             return None
169    
170      def getSegmentNormals(self, tag=None):
171         """
172         returns the list of the normals of the segments in fault ``tag``
173         :param tag: the tag of the fault
174         :type tag: ``float`` or ``str``
175         :return: the list of vectors normal to the segments. In the 2D case None is returned.
176         """
177         if tag==None: tag=self.NOTAG
178         if self.getDim()==3:
179             return self.__normals[tag]
180         else:
181             return None
182    
183      def getDepthVectors(self, tag=None):
184         """
185         returns the list of the depth vector at top vertices in fault ``tag``.
186         :param tag: the tag of the fault
187         :type tag: ``float`` or ``str``
188         :return: the list of segment depths. In the 2D case None is returned.
189         """
190         if tag==None: tag=self.NOTAG
191         if self.getDim()==3:
192             return self.__depth_vectors[tag]
193         else:
194             return None
195      def getDepths(self, tag=None):
196         """
197         returns the list of the depths of the segements in fault ``tag``.
198         :param tag: the tag of the fault
199         :type tag: ``float`` or ``str``
200         :return: the list of segment depths. In the 2D case None is returned.
201         """
202         if tag==None: tag=self.NOTAG
203         if self.getDim()==3:
204             return self.__depths[tag]
205         else:
206             return None
207    
208    def getW0Range(self,tag=None):    def getW0Range(self,tag=None):
209       """       """
# Line 105  class FaultSystem: Line 231  class FaultSystem:
231       return self.__offsets[tag]       return self.__offsets[tag]
232    
233    
   def getFaultSegments(self, tag=None):  
      """  
      returns the polygons used to describe fault tagged by ``tag``  
         
      :param tag: the tag of the fault  
      :type tag: ``float`` or ``str``  
      :return: the list of vertices defining the top and the list of the vertices defining the bottom of the fault indexed by ``tag``. The coordinates are `numpy.ndarray'.  
      """  
      if tag==None: tag=self.NOTAG  
      return self.__faults[tag]  
   
234    def getCenterOnSurface(self):    def getCenterOnSurface(self):
235        """        """
236        returns the center point of the fault system at the surface        returns the center point of the fault system at the surface
# Line 125  class FaultSystem: Line 240  class FaultSystem:
240          self.__center=numpy.zeros((3,),numpy.float)          self.__center=numpy.zeros((3,),numpy.float)
241          counter=0          counter=0
242          for t in self.getTags():          for t in self.getTags():
243              for s in self.getFaultSegments(t)[0]:              for s in self.getTopPolyline(t):
244                  self.__center[:2]+=s[:2]                  self.__center[:2]+=s[:2]
245                  counter+=1                  counter+=1
246          self.__center/=counter          self.__center/=counter
# Line 140  class FaultSystem: Line 255  class FaultSystem:
255            center=self.getCenterOnSurface()            center=self.getCenterOnSurface()
256            covariant=numpy.zeros((2,2))            covariant=numpy.zeros((2,2))
257            for t in self.getTags():            for t in self.getTags():
258                for s in self.getFaultSegments(t)[0]:                for s in self.getTopPolyline(t):
259                  covariant[0,0]+=(center[0]-s[0])**2                  covariant[0,0]+=(center[0]-s[0])**2
260                  covariant[0,1]+=(center[1]-s[1])*(center[0]-s[0])                  covariant[0,1]+=(center[1]-s[1])*(center[0]-s[0])
261                  covariant[1,1]+=(center[1]-s[1])**2                  covariant[1,1]+=(center[1]-s[1])**2
262                  covariant[1,0]+=(center[1]-s[1])*(center[0]-s[0])                  covariant[1,0]+=(center[1]-s[1])*(center[0]-s[0])
263            e, V=numpy.linalg.eigh(covariant)            e, V=numpy.linalg.eigh(covariant)
264            if e[0]>e[1]:            if e[0]>e[1]:
              d=V[:,1]  
           else:  
265               d=V[:,0]               d=V[:,0]
266              else:
267                 d=V[:,1]
268            if abs(d[0])>0.:            if abs(d[0])>0.:
269               self.__orientation=atan(d[1]/d[0])               self.__orientation=atan(d[1]/d[0])
270            else:            else:
# Line 168  class FaultSystem: Line 283  class FaultSystem:
283          mat=numpy.array([[cos(rot), -sin(rot)], [sin(rot), cos(rot)] ])          mat=numpy.array([[cos(rot), -sin(rot)], [sin(rot), cos(rot)] ])
284       else:       else:
285          mat=numpy.array([[cos(rot), -sin(rot),0.], [sin(rot), cos(rot),0.], [0.,0.,1.] ])          mat=numpy.array([[cos(rot), -sin(rot),0.], [sin(rot), cos(rot),0.], [0.,0.,1.] ])
286    
287       for t in self.getTags():       for t in self.getTags():
288             top=[]           strikes=[ s+ rot for s in self.getStrikes(t) ]
289             for s in self.getFaultSegments(t)[0]: top.append(numpy.dot(mat,(s+shift[:self.getDim()])))           V0=self.getStart(t)
290             if self.getDim() == 3:  
291                bottom=[]           self.addFault(strikes = [ s+ rot for s in self.getStrikes(t) ], \
292                for s in self.getFaultSegments(t)[1]: bottom.append(numpy.dot(mat,(s+shift[:self.getDim()])))                         ls = self.getLengths(t), \
293             else:                         V0=numpy.dot(mat,self.getStart(t)+shift), \
294                bottom=None                         tag =t, \
295             self.addFault(top=top, tag=t, bottom=bottom, w0_offsets=self.getW0Offsets(t), w1_max=-self.getW1Range(t)[0])                         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, top, tag=None, bottom=None, w0_offsets=None, w1_max=None):    def addFault(self, strikes, ls, V0=[0.,0.,0.],tag=None, dips=None, depths= None, w0_offsets=None, w1_max=None):
301       """       """
302       adds a new fault to the fault system. The fault is named by ``tag``.       adds a new fault to the fault system. The fault is named by ``tag``.
303    
304         The fault is defined by a starting point V0 and a list of ``strikes`` and length ``ls``. The strikes and the length
305         is used to define a polyline with points ``V[i]`` such that
306    
307       a fault is subdivided into segments where segment ``i`` is represented by the four corner points       - ``V[0]=V0``
308       ``v0``=``top[i]``, ``v1``=``bottom[i]``, ``v2``=``bottom[i+1]`` and ``v3``==``top[i+1]``. The segment is parametrized       - ``V[i]=V[i]+ls[i]*array(cos(strikes[i]),sin(strikes[i]),0)''
309       by ``w0`` and ``w1`` with ``w0_offsets[i]<=w0<=w0_offsets[i+1]`` and ``-w1_max<=w1<=0``. Moreover  
310         So ``strikes`` defines the angle between the direction of the fault segment and the x0 axis. In 3D ''ls[i]``==0 is allowed.
311    
312         In case of a 3D model a fault plane is defined through a dip ``dips`` and depth ``depths``.
313         From the dip and the depth the polyline ``bottom`` of the bottom of the fault is computed.
314    
315    
316         Each segment in the fault is decribed by the for vertices ``v0=top[i]``, ``v1``==``top[i+1]``, ``v2=bottom[i]`` and ``v3=bottom[i+1]``
317         The segment is parametrized by ``w0`` and ``w1`` with ``w0_offsets[i]<=w0<=w0_offsets[i+1]`` and ``-w1_max<=w1<=0``. Moreover
318        
319       - ``(w0,w1)=(w0_offsets[i]  ,       0)->v0``       - ``(w0,w1)=(w0_offsets[i]  ,       0)->v0``
320       - ``(w0,w1)=(w0_offsets[i]  , -w1_max)->v1``       - ``(w0,w1)=(w0_offsets[i+1],       0)->v1``
321       - ``(w0,w1)=(w0_offsets[i+1], -w1_max)->v2``       - ``(w0,w1)=(w0_offsets[i]  , -w1_max)->v2``
322       - ``(w0,w1)=(w0_offsets[i+1],       0)->v3``       - ``(w0,w1)=(w0_offsets[i+1], -w1_max)->v3``
323    
324       If no ``w0_offsets`` is given,       If no ``w0_offsets`` is given,
325        
326       - ``w0_offsets[0]=0``       - ``w0_offsets[0]=0``
327       - ``w0_offsets[i]=w0_offsets[i-1]+length(v0-v4)``       - ``w0_offsets[i]=w0_offsets[i-1]+L[i]``
328    
329       If no ``w1_max`` is given, the average fault depth (=length(v0-v1)) is used.       where ``L[i]`` is the length of the segments on the top in 2D and in the middle of the segment in 3D.
330    
331       :param top: list of coordinates. top[i] and top[i+1] are the top corners of the i-th segment of the fault       If no ``w1_max`` is given, the average fault depth is used.
332       :type top: ``list`` of objects which can be converted to `numpy.ndarray` vector of length 3 (or 2 in the 2D case)  
333    
334         :param strikes: list of strikes. This is the angle of the fault segment direction with x0 axis. Right hand rule applies.
335         :type strikes: ``list`` of ``float``
336         :param ls: list of fault lengths. In the case of a 3D fault a segment may have length 0.
337         :type ls: ``list`` of ``float``
338         :param V0: start point of the fault
339         :type V0: ``list`` or `numpy.ndarray' with 2 or 3 components. ``V0[2]`` must be zero.
340       :param tag: the tag of the fault. If fault ``tag`` already exists it is overwritten.       :param tag: the tag of the fault. If fault ``tag`` already exists it is overwritten.
341       :type tag: ``float`` or ``str``       :type tag: ``float`` or ``str``
342       :param bottom: list of coordinates. bottom[i] and bottom[i+1] are the bottom corners of the i-th segment of the fault       :param dips: list of dip angles. Right hand rule around strike direction applies.
343       :type bottom: ``list`` of objects which can be converted to `numpy.ndarray` vector of length 3 (or 2 in the 2D case)       :type dips: ``list`` of ``float``
344         :param depths: list of segment depth. Value mut be positive in the 3D case.
345         :type depths: ``list`` of ``float``
346       :param w0_offsets: ``w0_offsets[i]`` defines the offset of the segment ``i`` in the fault to be used in the parametrization of the fault. If not present the cumulative length of the fault segments is used.       :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``       :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.       :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``       :type w1_max: ``float``
350       :note: In the three dimensional case the lists ``bottom`` and ``top`` must have the same length.       :note: In the three dimensional case the lists ``dip`` and ``top`` must have the same length.
351       """       """
352       if tag==None:       if tag==None:
353           tag=self.NOTAG           tag=self.NOTAG
354       else:       else:
355           if self.NOTAG in self.getTags():           if self.NOTAG in self.getTags():
356                raise ValueError,'Attempt to add a fault with no tag to a set of existing faults'                raise ValueError,'Attempt to add a fault with no tag to a set of existing faults'
357       if len(top)<2:       if not isinstance(strikes, list): strikes=[strikes, ]
358           raise Value,"at least two vertices must be given."       n_segs=len(strikes)
359       if self.getDim()==2 and not bottom==None:       if not isinstance(ls, list): ls=[ ls for i in xrange(n_segs) ]
360             raise ValueError,'Spatial dimension two does not support bottom polygon for faults.'       if not n_segs==len(ls):
361       if not bottom == None:           raise ValueError,"number of strike direction and length must match."
362          if len(top) != len(bottom):       if len(V0)>2:
363             raise ValueError,'length of top and bottom polygon must be identical.'            if abs(V0[2])>0: raise Value,"start point needs to be surface (3rd component ==0)"
364         if self.getDim()==2 and not  (dips==None and depths == None) :
365               raise ValueError,'Spatial dimension two does not support dip and depth for faults.'
366         if not dips == None:
367            if not isinstance(dips, list): dips=[dips for i in xrange(n_segs) ]
368            if n_segs != len(dips):
369               raise ValueError,'length of dips must be one less then the length of top.'
370         if not depths == None:
371            if not isinstance(depths, list): depths=[depths for i in xrange(n_segs+1) ]
372            if n_segs+1 != len(depths):
373               raise ValueError,'length of depths must be one less then the length of top.'
374       if w0_offsets != None:       if w0_offsets != None:
375         if len(w0_offsets) != len(top):         if len(w0_offsets) != n_segs+1:
376            raise ValueError,'expected length of w0_offsets is %s'%(len(top))            raise ValueError,'expected length of w0_offsets is %s'%(n_segs)
377       self.__center=None       self.__center=None
378       self.__orientation = None       self.__orientation = None
      #  
      #   translate to numpy  
379       #       #
380       if self.getDim()==2:       #  in the 2D case we don't allow zero length:
381          self.__faults[tag]=([ numpy.array(t[:2],numpy.double) for t in top ],)       #
382         if self.getDim() == 2:
383            for l in ls:
384                if l<=0: raise ValueError,"length must be positive"
385       else:       else:
386          self.__faults[tag]=( [ numpy.array(t[:3],numpy.double) for t in top ], [ numpy.array(b[:3],numpy.double) for b in bottom ])          for l in ls:
387       faults=self.__faults[tag]              if l<0: raise ValueError,"length must be non-negative"
388       len_faults=len(faults[0])          for i in xrange(n_segs+1):
389               if depths[i]<0: raise ValueError,"negative depth."
390       #       #
391       #    calculate the depth       #   translate start point to numarray
392       #       #
393       if self.getDim()==2:       V0= numpy.array(V0[:self.getDim()],numpy.double)
         self.__depth[tag]=0.  
      else:  
         self.__depth[tag]=float(sum([numpy.linalg.norm(faults[1][i]-faults[0][i]) for i in xrange(len_faults)])/len_faults)  
      if w1_max==None or self.getDim()==2: w1_max=self.__depth[tag]  
      self.__length[tag]=float(sum([numpy.linalg.norm(faults[0][i]-faults[0][i-1]) for i in xrange(1,len_faults)]))  
      self.__w1_max[tag]=w1_max  
394       #       #
395       #       #  set strike vectors:
396       #       #
397       if w0_offsets!=None:       strike_vectors=[]
398          self.__offsets[tag]=w0_offsets       top_polyline=[V0]  
399         total_length=0
400         for i in xrange(n_segs):
401             v=numpy.zeros((self.getDim(),))
402             v[0]=cos(strikes[i])
403             v[1]=sin(strikes[i])
404             strike_vectors.append(v)
405             top_polyline.append(top_polyline[-1]+ls[i]*v)
406             total_length+=ls[i]
407         #
408         #    normal and depth direction
409         #
410         if self.getDim()==3:
411            normals=[]
412            for i in xrange(n_segs):
413               normals.append(numpy.array([sin(dips[i])*strike_vectors[i][1], -sin(dips[i])*strike_vectors[i][0], cos(dips[i])]) )
414      
415            d=numpy.cross(strike_vectors[0],normals[0])
416            depth_vectors=[depths[0]*d/numpy.linalg.norm(d) ]
417            for i in xrange(1,n_segs):
418                d=-numpy.cross(normals[i-1],normals[i])
419                d_l=numpy.linalg.norm(d)
420                if d_l<=0:
421                     d=numpy.cross(strike_vectors[i],normals[i])
422                     depth_vectors.append([depths[d]*d/numpy.linalg.norm(d)])
423                else:
424                     d=depths[i]*d/d_l
425                     for L in [ strike_vectors[i], strike_vectors[i-1]]:
426                        if numpy.linalg.norm(numpy.cross(L,d)) <= self.MIN_DEPTH_ANGLE * numpy.linalg.norm(L) * depths[i]:
427                             raise ValueError,"%s-th depth vector %s too flat."%(i, d)
428                     depth_vectors.append(d)
429            d=numpy.cross(strike_vectors[n_segs-1],normals[n_segs-1])
430            depth_vectors.append(depths[n_segs]*d/numpy.linalg.norm(d))
431            bottom_polyline=[ top_polyline[i]+depth_vectors[i] for i in xrange(n_segs+1) ]
432         #
433         #   calculate offsets if required:
434         #
435         if w0_offsets==None:
436            w0_offsets=[0.]
437            for  i in xrange(n_segs):
438                if self.getDim()==3:
439                   w0_offsets.append(w0_offsets[-1]+(float(numpy.linalg.norm(bottom_polyline[i+1]-bottom_polyline[i]))+ls[i])/2.)
440                else:
441                   w0_offsets.append(w0_offsets[-1]+ls[i])
442         w0_max=max(w0_offsets)
443         if self.getDim()==3:
444            self.__normals[tag]=normals
445            self.__depth_vectors[tag]=depth_vectors
446            self.__depths[tag]=depths
447            self.__dips[tag]=dips
448            self.__bottom[tag]=bottom_polyline
449         self.__ls[tag]=ls
450         self.__strikes[tag]=strikes
451         self.__strike_vectors[tag]=strike_vectors
452         self.__top[tag]=top_polyline
453         self.__total_length[tag]=total_length
454         self.__offsets[tag]=w0_offsets
455    
456         if self.getDim()==2:
457            self.__medDepth[tag]=0.
458       else:       else:
459          self.__offsets[tag]=[0.]          self.__medDepth[tag]=sum([ numpy.linalg.norm(v) for v in depth_vectors])/len(depth_vectors)
460          for  i in xrange(1,len_faults):       if w1_max==None or self.getDim()==2: w1_max=self.__medDepth[tag]
             self.__offsets[tag].append(self.__offsets[tag][-1]+float(numpy.linalg.norm(faults[0][i]-faults[0][i-1])))  
      w0_max=max(self.__offsets[tag])  
461       self.__w0_max[tag]=w0_max       self.__w0_max[tag]=w0_max
462         self.__w1_max[tag]=w1_max
463    
464    def getMaxValue(self,f, tol=sqrt(EPSILON)):    def getMaxValue(self,f, tol=sqrt(EPSILON)):
465       """       """
# Line 301  class FaultSystem: Line 504  class FaultSystem:
504    
505      :param x: location(s)      :param x: location(s)
506      :type x: `escript.Data` object or `numpy.ndarray`      :type x: `escript.Data` object or `numpy.ndarray`
507      :param tag: the tage of the fault      :param tag: the tag of the fault
508      :param tol: relative tolerance to check if location is on fault.      :param tol: relative tolerance to check if location is on fault.
509      :type tol: ``float``      :type tol: ``float``
510      :param outsider: value used for parametrization values outside the fault. If not present an appropriate value is choosen.      :param outsider: value used for parametrization values outside the fault. If not present an appropriate value is choosen.
# Line 310  class FaultSystem: Line 513  class FaultSystem:
513      :rtype: `escript.Data` object or `numpy.ndarray`      :rtype: `escript.Data` object or `numpy.ndarray`
514      """      """
515      offsets=self.getW0Offsets(tag)      offsets=self.getW0Offsets(tag)
516        w1_range=self.getW1Range(tag)
517      w0_range=self.getW0Range(tag)[1]-self.getW0Range(tag)[0]      w0_range=self.getW0Range(tag)[1]-self.getW0Range(tag)[0]
518      if outsider == None:      if outsider == None:
519         outsider=min(self.getW0Range(tag)[0],self.getW0Range(tag)[1])-abs(w0_range)/sqrt(EPSILON)         outsider=min(self.getW0Range(tag)[0],self.getW0Range(tag)[1])-abs(w0_range)/sqrt(EPSILON)
520                    
521      p=x[0]*0 + outsider      if isinstance(x,list): x=numpy.array(x, numpy.double)
522      updated=x[0]*0      updated=x[0]*0
523    
524      if self.getDim()==2:      if self.getDim()==2:
525          #          #
526          #          #
527          s=self.getFaultSegments(tag)[0]          p=x[0]*0 + outsider
528          for i in xrange(1,len(s)):          top=self.getTopPolyline(tag)
529             d=s[i]-s[i-1]          for i in xrange(1,len(top)):
530             h=x-s[i-1]             d=top[i]-top[i-1]
531               h=x-top[i-1]
532             h_l=length(h)             h_l=length(h)
533             d_l=length(d)             d_l=length(d)
534             if not d_l>0:             s=inner(h,d)/d_l**2
535                raise ValueError,"segement %s in fault %s has zero length."%(i,tag)             s=s*whereNonPositive(s-1.-tol)*whereNonNegative(s+tol)
536             t=inner(h,d)/d_l**2             m=whereNonPositive(length(h-s*d)-tol*maximum(h_l,d_l))*(1.-updated)
537             t=t*whereNonPositive(t-1.)*whereNonNegative(t)             p=(1.-m)*p+m*(offsets[i-1]+(offsets[i]-offsets[i-1])*s)
            r=length(h-t*d)/maximum(h_l,d_l)  
            m=whereNonPositive(r-tol)*(1.-updated)  
            p=(1.-m)*p+m*(offsets[i-1]+(offsets[i]-offsets[i-1])*t)  
538             updated=wherePositive(updated+m)             updated=wherePositive(updated+m)
539      else:      else:
540         raise ValueError,"3D is not supported yet."          p=x[:2]*0 + outsider
541            top=self.getTopPolyline(tag)
542            bottom=self.getBottomPolyline(tag)
543            n=self.getSegmentNormals(tag)
544            for i in xrange(len(top)-1):
545                h=x-top[i]
546                R=top[i+1]-top[i]
547                r=bottom[i+1]-bottom[i]
548                D0=bottom[i]-top[i]
549                D1=bottom[i+1]-top[i+1]
550                s_upper=matrix_mult(numpy.linalg.pinv(numpy.vstack((R,D1)).T),h)
551                s_lower=matrix_mult(numpy.linalg.pinv(numpy.vstack((r,D0)).T),h)
552                m_ul=wherePositive(s_upper[0]-s_upper[1])
553                s=s_upper*m_ul+s_lower*(1-m_ul)
554                s0=s[0]
555                s1=s[1]
556                m=whereNonNegative(s0+tol)*whereNonPositive(s0-1.-tol)*whereNonNegative(s1+tol)*whereNonPositive(s1-1.-tol)
557                s0=s0*m
558                s1=s1*m
559                atol=tol*maximum(length(h),length(top[i]-bottom[i+1]))
560                m=whereNonPositive(length(h-s0*R-s1*D1)*m_ul+(1-m_ul)*length(h-s0*r-s1*D0)-atol)
561                p[0]=(1.-m)*p[0]+m*(offsets[i]+(offsets[i+1]-offsets[i])*s0)
562                p[1]=(1.-m)*p[1]+m*(w1_range[1]+(w1_range[0]-w1_range[1])*s1)
563                updated=wherePositive(updated+m)
564            
565      return p, updated      return p, updated
566    
# Line 345  class FaultSystem: Line 570  class FaultSystem:
570    
571      :param x: location(s)      :param x: location(s)
572      :type x: `escript.Data` object or `numpy.ndarray`      :type x: `escript.Data` object or `numpy.ndarray`
573      :param tag: the tage of the fault      :param tag: the tag of the fault
574      :return: the side of ``x`` (positive means to the right of the fault, negative to the left) and the distance to the fault. Note that a value zero for the side means that that the side is undefined.      :return: the side of ``x`` (positive means to the right of the fault, negative to the left) and the distance to the fault. Note that a value zero for the side means that that the side is undefined.
575      """      """
576      d=None      d=None
577      side=None      side=None
578      if self.getDim()==2:      if self.getDim()==2:
579          mat=numpy.array([[0., 1.], [-1., 0.] ])          mat=numpy.array([[0., 1.], [-1., 0.] ])
580          s=self.getFaultSegments(tag)[0]          s=self.getTopPolyline(tag)
581          for i in xrange(1,len(s)):          for i in xrange(1,len(s)):
582             q=(s[i]-s[i-1])             q=(s[i]-s[i-1])
583             h=x-s[i-1]             h=x-s[i-1]
584             q_l=length(q)             q_l=length(q)
            if not q_l>0:  
               raise ValueError,"segement %s in fault %s has zero length."%(i,tag)  
585             qt=matrixmult(mat,q)   # orthogonal direction             qt=matrixmult(mat,q)   # orthogonal direction
586             t=inner(q,h)/q_l**2             t=inner(q,h)/q_l**2
587             t=maximum(minimum(t,1,),0.)             t=maximum(minimum(t,1,),0.)
# Line 374  class FaultSystem: Line 597  class FaultSystem:
597                 d=dist*(1-m)+d*m                 d=dist*(1-m)+d*m
598                 side=lside*(1-m2)+side*m2                 side=lside*(1-m2)+side*m2
599      else:      else:
600         raise ValueError,"3D is not supported yet."          ns=self.getSegmentNormals(tag)
601      return side, d          top=self.getTopPolyline(tag)
602            bottom=self.getBottomPolyline(tag)
603            for i in xrange(len(top)-1):
604                h=x-top[i]
605                R=top[i+1]-top[i]
606                r=bottom[i+1]-bottom[i]
607                D0=bottom[i]-top[i]
608                D1=bottom[i+1]-top[i+1]
609                s_upper=matrix_mult(numpy.linalg.pinv(numpy.vstack((R,D1)).T),h)
610                s_lower=matrix_mult(numpy.linalg.pinv(numpy.vstack((r,D0)).T),h)
611                m_ul=wherePositive(s_upper[0]-s_upper[1])
612                s=s_upper*m_ul+s_lower*(1-m_ul)
613                s=maximum(minimum(s,1.),0)
614                p=h-s[0]*(R*m_ul+(1-m_ul)*r)-s[1]*(D1*m_ul+(1-m_ul)*D0)
615                dist=length(p)
616                lside=sign(inner(p,ns[i]))
617                if d == None:
618                   d=dist
619                   side=lside
620                else:
621                   m=whereNegative(d-dist)
622                   m2=wherePositive(whereZero(abs(lside))+m)
623                   d=dist*(1-m)+d*m
624                   side=lside*(1-m2)+side*m2
625    
626  def patchMap(v,v1,v2,v3,v4,x_offset,x_length,depth):      return side, d
     """  
     maps a patch with corners v1,v2,v3,v4 onto  
     the rectangle [x_offset,0] x [x_offset+x_length,-depth]  
     """  
     #  
     #   inspect the mapping v-> w=(w0,w1) with 0<=s<=1 0<=t<=1  
     #   then we apply x_offset+s*x_length, -depth*t to map to the final domain  
     #  
     #   we need to find s and t such that  
     #  
     #   v=v1*(1-t)*(1-s)+v2*t*(1-s)+v3*t*s+v4*(1-t)*s  
     #    =v1*(1-t-s+t*s)+v2*(t-t*s)+v3*t*s+v4*(s-t*s)  
     #    =v1+t*(v2-v1)+s*(v4-v1)+s*t*(v1-v2+v3-v4)  
     #    
     #   with d4=(v2-v1)-dot(v2-v1,v4-v1)/dot(v4-v1,v4-v1)*(v4-v1) => dot(d4,v2-v1)=0  
     #   with d2=(v4-v1)-dot(v2-v1,v4-v1)/dot(v2-v1,v2-v1)*(v2-v1) => dot(d2,v4-v1)=0  
     #  
     #   dot(v-v1,d4) = s*dot(v4-v1,d4)+s*t*dot(v3-v4,d4)   (a)  
     #   dot(v-v1,d2) = t*dot(v2-v1,d2)+s*t*dot(v3-v2,d2)   (b)  
     #  
     #   in addition we set n=cross(d2,d4) => dot(v2-v1,n)=dot(v4-v1,n)=0  
     #  
     #   dot(v-v1,n)=s*t*dot(v3-v1,n)  
     #  
     #   => s*t=dot(v-v1,n)/dot(v3-v1,n) if dot(v3-v1,n)!=0  
     #      otherwise on can set s*t=0 to solve (a) and (b)  
     #      
     h=v-v1  
     h2=(v2-v1)  
     h3=(v3-v1)  
     h4=(v4-v1)  
     d4=h4-numpy.dot(h2,h4)/numpy.dot(h2,h2)*h2  
     d2=h2-numpy.dot(h2,h4)/numpy.dot(h4,h4)*h4  
     n=numpy.cross(d4,d2)  
     h31n=numpy.dot(n,v3-v1)  
     if abs(h31n)>0:  
         st=n/h31n  
     else:  
         st=0*n  
     s=numpy.dot(h,d4-st*numpy.dot(v3-v4,d4))/numpy.dot(v4-v1,d4)  
     t=numpy.dot(h,d2-st*numpy.dot(v3-v2,d2))/numpy.dot(v2-v1,d2)  
     return x_offset+s*x_length, -depth*t  
627    

Legend:
Removed from v.2662  
changed lines
  Added in v.2663

  ViewVC Help
Powered by ViewVC 1.1.26