/[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 2653 by gross, Mon Sep 7 03:33:55 2009 UTC revision 2654 by gross, Tue Sep 8 07:11:12 2009 UTC
# Line 339  class FaultSystem: Line 339  class FaultSystem:
339            
340      return p, updated      return p, updated
341    
342  def __getSignedDistance(self,x,tag=None):    def getSideAndDistance(self,x,tag=None):
343      """      """
344      returns the signed distance at ``x`` from the fault ``tag`` where the first and the last segments are extended to infinity.      returns the side and the distance at ``x`` from the fault ``tag``.
345    
346      :param x: location(s)      :param x: location(s)
347      :type x: `escript.Data` object or `numpy.ndarray`      :type x: `escript.Data` object or `numpy.ndarray`
348      :param tag: the tage of the fault      :param tag: the tage of the fault
349        :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.
350      """      """
351      p=x[0]*0      d=None
352        side=None
353      if self.getDim()==2:      if self.getDim()==2:
354          mat=numpy.array([[0., -1.], [1., 0.] ])          mat=numpy.array([[0., 1.], [-1., 0.] ])
         #  
         #  
355          s=self.getFaultSegments(tag)[0]          s=self.getFaultSegments(tag)[0]
356          for i in xrange(1,len(s)):          for i in xrange(1,len(s)):
357             d=(s[i]-s[i-1])             q=(s[i]-s[i-1])
358             h=x-s[i-1]             h=x-s[i-1]
359             d_l=length(d)             q_l=length(q)
360             if not d_l>0:             if not q_l>0:
361                raise ValueError,"segement %s in fault %s has zero length."%(i,tag)                raise ValueError,"segement %s in fault %s has zero length."%(i,tag)
362             dt=numpy.dot(mat,d)/d_l             qt=matrixmult(mat,q)   # orthogonal direction
363             a=numpy.dot(h,dt)             t=inner(q,h)/q_l**2
364             p=maximum(a,p)             t=maximum(minimum(t,1,),0.)
365               p=h-t*q
366               dist=length(p)
367               lside=sign(inner(p,qt))
368               if d == None:
369                   d=dist
370                   side=lside
371               else:
372                   m=whereNegative(d-dist)
373                   m2=wherePositive(whereZero(abs(lside))+m)
374                   d=dist*(1-m)+d*m
375                   side=lside*(1-m2)+side*m2
376      else:      else:
377         raise ValueError,"3D is not supported yet."         raise ValueError,"3D is not supported yet."
378      return p      return side, d
   
379    
380  def patchMap(v,v1,v2,v3,v4,x_offset,x_length,depth):  def patchMap(v,v1,v2,v3,v4,x_offset,x_length,depth):
381      """      """

Legend:
Removed from v.2653  
changed lines
  Added in v.2654

  ViewVC Help
Powered by ViewVC 1.1.26