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

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

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

revision 1716 by gross, Thu Aug 21 05:03:49 2008 UTC revision 1788 by gross, Mon Sep 15 02:20:20 2008 UTC
# Line 6  import math Line 6  import math
6    
7    
8    
9  class LevelSet(object):  from esys.escript import *
10       def __init__(self,phi,reinitialization_steps_max=10,reinitialization_each=2,relative_smoothing_width=2.):  from esys.finley import finley
11    from esys.escript.linearPDEs import LinearPDE
12    from esys.escript.pdetools import Projector
13    import sys
14    
15    
16    class LevelSet:
17      """
18      The level set method tracking an interface defined by the zero contour of the level set function phi.
19      it is assumed the phi(x)<0 defines the volume of interest.
20    
21      """
22      def __init__(self,domain,phi,reinit_max=10,reinit_each=2,smooth=2.):
23        """
24        set up the level set method
25    
26        @param domain: the domain where the level set is used
27        @param phi: the initial level set function
28        @param reinit_max: maximum number of reinitalization steps
29        @param reinit_each: phi is reinitialized every reinit_each step
30        @param smooth: smoothing width
31        """
32        self.__domain = domain
33        self.__phi = phi
34        self.__reinit_max = reinit_max
35        self.__reinit_each = reinit_each
36        self.__PDE = LinearPDE(domain)
37        self.__PDE.setReducedOrderOn()
38        self.__PDE.setValue(D=1.0)
39        self.__PDE.setSolverMethod(solver=LinearPDE.PCG)
40        self.__reinitPDE = LinearPDE(domain, numEquations=1)
41        self.__reinitPDE.setReducedOrderOn()
42        self.__reinitPDE.setValue(D=1.0)
43        self.__reinitPDE.setSolverMethod(solver=LinearPDE.LUMPING)
44        self.__h = inf(domain.getSize())
45        self.__smooth = smooth
46        self.__n_step=0
47      
48      def __advect(self, velocity, dt):
49        """
50        advects the level set function in the presense of a velocity field.
51    
52        This implementation uses the 2-step Taylor-Galerkin method
53        @param velocity: velocity field
54        @param dt: dime increment
55        @return: the advected level set function
56        """
57        Y = self.__phi-dt/2.*inner(velocity,grad(self.__phi))
58        self.__PDE.setValue(Y=Y)    
59        phi_half = self.__PDE.getSolution()
60        Y = self.__phi-dt*inner(velocity,grad(phi_half))
61        self.__PDE.setValue(Y=Y)    
62        phi = self.__PDE.getSolution()
63        print "LevelSet: Advection step done"
64        return phi
65    
66      def __reinitialise(self):
67        """
68        reinializes the level set
69    
70        It solves the
71    
72        @return: reinitalized level set
73        """
74        phi=self.__phi
75        s = sign(phi.interpolate(Function(self.__domain)))
76        g=grad(phi)
77        w = s*g/(length(g)+1e-10)
78        dtau = 0.5*self.__h
79        iter =0
80        mask = whereNegative(abs(phi)-1.*self.__h)
81        self.__reinitPDE.setValue(q=mask, r=phi)
82        g=grad(phi)
83        while (iter<=self.__reinit_max):
84          Y = phi+(dtau/2.)*(s-inner(w,g))
85          self.__reinitPDE.setValue(Y = Y)
86          phi_half = self.__reinitPDE.getSolution()
87          Y = phi+dtau*(s-inner(w,grad(phi_half)))
88          self.__reinitPDE.setValue(Y = Y)
89          phi = self.__reinitPDE.getSolution()
90          g=grad(phi)
91          error = Lsup(length(g)*whereNegative(abs(phi.interpolate(g.getFunctionSpace()))-3.0*self.__h))
92          print "LevelSet:reinitialization: iteration :", iter, " error:", error
93          iter +=1
94        return phi
95    
96      def getTimeStepSize(self,velocity):
97           """
98           returns a new dt for a given velocity using the courant coundition
99    
100           @param velocity: velocity field
101           """
102           self.__velocity=velocity
103           dt=0.5*self.__h/sup(length(velocity))
104           return dt
105      
106      def update(self,dt):
107          """
108          sets a new velocity and updates the level set fuction
109    
110          @param dt: time step forward
111          """
112          self.__phi=self.__advect(self.__velocity, dt)
113          self.__n_step+=1
114          if self.__n_step%self.__reinit_each ==0: self.__phi = self.__reinitialise()
115          return self.__phi
116    
117      def update_phi(self, velocity, dt):
118          """
119          updates phi under the presense of a velocity field
120    
121          If dt is small this call is equivalent to call
122      
123          dt=LevelSet.getTimeStepSize(velocity)
124          phi=LevelSet.update(dt)
125                                      
126          otherwise substepping is used.
127          @param velocity: velocity field
128          @param dt: time step forward
129          """
130          dt2=self.getTimeStepSize(velocity)
131          n=math.ceil(dt/dt2)
132          dt_new=dt/n
133          for i in range(n):
134               phi=self.update(dt_new)
135               t+=dt_new
136          return phi
137    
138    
139      def getVolume(self):
140        """
141        return the volume of the phi(x)<0 region
142        """
143        return integrate(whereNegative(self.__phi.interpolate(Function(self.__domain))))
144    
145      def getSurface(self,rel_width_factor=0.5):
146        """
147        return a mask for phi(x)=1 region
148      
149        @param rel_width_factor: relative wideth of region around zero contour.
150        """
151        return whereNegative(abs(self.__phi)-rel_width_factor*self.__h)
152        
153      def getH(self):
154         """
155         returns mesh size
156         """
157         return self.__h
158    
159      def getDomain(self):
160         """
161         returns domain
162         """
163         return self.__domain
164    
165      def getLevelSetFunction(self):
166          """
167          returns the level set function
168          """
169          return self.__phi
170    
171      def update_parameter_sharp(self, param_neg=-1, param_pos=1, phi=None):
172        """
173        creates a function whith param_neg where phi<0 and param_pos where phi>0 (no smoothing)
174    
175        @param param_neg: value of parameter on the negative side (phi<0)
176        @param param_pos: value of parameter on the positve side (phi>0)
177        @param phi: level set funtion to be used. if not present the current level set is used.
178        """
179        mask_neg = whereNegative(self.__phi)
180        mask_pos = whereNonNegative(self.__phi)
181        param = param_pos*mask_pos + param_neg*mask_neg
182        return param
183    
184      def update_parameter(self, param_neg=-1, param_pos=1, phi=None, smoothing_width=None):
185        """
186        creates a smoothed function whith param_neg where phi<0 and param_pos where phi>0 which is smoothed over a length
187        smoothing_width accross the interface
188    
189        @param smoothing_width: width of the smoothing zone relative to mesh size. If not present the initial value of C{smooth} is used.
190        """
191        if smoothing_width==None: smoothing_width = self.__smooth
192        if phi==None: phi=self.__phi
193        s=self.__makeInterface(phi,smoothing_width)
194        return ((param_pos-param_neg)*s+param_pos+param_neg)/2
195    
196      def __makeInterface(self,phi,smoothing_width):
197          """
198          creates a smooth interface from -1 to 1 over the length 2*h*smoothing_width where -1 is used where the level set is negative
199          and 1 where the level set is 1
200          """
201          s=smoothing_width*self.__h
202          phi_on_h=interpolate(phi,Function(self.__domain))        
203          mask_neg = whereNonNegative(-s-phi_on_h)
204          mask_pos = whereNonNegative(phi_on_h-s)
205          mask_interface = 1.-mask_neg-mask_pos
206          interface=phi_on_h/s
207          return - mask_neg + mask_pos + mask_interface * interface
208    
209      def makeCharacteristicFunction(self, contour=0, phi=None, positiveSide=True, smoothing_width=None):
210          """
211          makes a smooth charateristic function of the region phi(x)>contour if positiveSide and phi(x)<contour otherwise.
212    
213          @param phi: level set funtion to be used. if not present the current level set is used.
214          @param smoothing_width: width of the smoothing zone relative to mesh size. If not present the initial value of C{smooth} is used.
215          """
216          if phi==None: phi=self.__phi
217          if smoothing_width == None: smoothing_width=self.__smooth
218          s=self.__makeInterface(phi=phi-contour,smoothing_width=smoothing_width)
219          if positiveSide:
220              return (1+s)/2
221          else:
222              return (1-s)/2
223    
224      def setTolerance(self,tolerance=1e-3):
225        self.__PDE.setTolerance(tolerance)
226        self.__reinitPDE.setTolerance(tolerance)
227    
228    
229    class LevelSet2(object):
230         def __init__(self,phi,reinit_max=10,reinit_each=2,smooth=2.):
231           """           """
232           initialize model           initialize model
233           """           """
# Line 19  class LevelSet(object): Line 239  class LevelSet(object):
239               diam+=(inf(xi)-sup(xi))**2               diam+=(inf(xi)-sup(xi))**2
240           self.__diam=sqrt(diam)           self.__diam=sqrt(diam)
241           self.__h = sup(Function(self.__domain).getSize())           self.__h = sup(Function(self.__domain).getSize())
242           self.__reinitialization_each=reinitialization_each           self.__reinit_each=reinit_each
243           self.__reinitialization_steps_max = reinitialization_steps_max           self.__reinit_max = reinit_max
244           self.__relative_smoothing_width = relative_smoothing_width           self.__smooth = smooth
245           self.__phi = phi           self.__phi = phi
246           self.__update_count=0           self.__update_count=0
247           self.velocity = None           self.velocity = None
# Line 63  class LevelSet(object): Line 283  class LevelSet(object):
283           # self.__fc.setValue(q=q)           # self.__fc.setValue(q=q)
284           self.__reinitPde.setValue(q=q)           self.__reinitPde.setValue(q=q)
285    
      def getH(self):  
          return self.__h  
      def getDomain(self):  
          return self.__domain  
   
      def getTimeStepSize(self,velocity):  
          """  
          returns a new dt for a given velocity using the courant coundition  
          """  
          self.velocity=velocity  
          if self.__FC:  
             self.__fc.setValue(C=-interpolate(velocity,Function(self.__domain)))  
             dt=self.__fc.getSafeTimeStepSize()  
          else:  
             dt=0.5*self.__h/sup(length(velocity))  
          return dt  
   
      def getLevelSetFunction(self):  
          """  
          returns the level set function  
          """  
          return self.__phi  
286    
287       def __updateInterface(self):       def __updateInterface(self):
288           self.__smoothed_char=self.__makeInterface(self.__phi)           self.__smoothed_char=self.__makeInterface(self.__phi)
289    
290       def __makeInterface(self,phi,smoothing_width=1.):  
          """  
          creates a very smooth interface from -1 to 1 over the length 2*h*smoothing_width where -1 is used where the level set is negative  
          and 1 where the level set is 1  
          """  
          s=smoothing_width*self.__h  
          phi_on_h=interpolate(phi,Function(self.__domain))          
          mask_neg = whereNonNegative(-s-phi_on_h)  
          mask_pos = whereNonNegative(phi_on_h-s)  
          mask_interface = 1.-mask_neg-mask_pos  
          # interface=1.-(phi_on_h-s)**2/(2.*s**3)*(phi_on_h+2*s)  # function f with f(s)=1, f'(s)=0, f(-s)=-1, f'(-s)=0, f(0)=0, f'(0)=  
          interface=phi_on_h/s  
          return - mask_neg + mask_pos + mask_interface * interface  
   
      def makeCharacteristicFunction(self, contour=0, positiveSide=True, smoothing_width=1.):  
          return self.makeCharacteristicFunctionFromExternalLevelSetFunction(self.__phi,contour,positiveSide,smoothing_width)  
   
      def makeCharacteristicFunctionFromExternalLevelSetFunction(self,phi,contour=0, positiveSide=True, smoothing_width=1.):  
          s=self.__makeInterface(phi-contour,smoothing_width)  
          if positiveSide:  
             return (1+s)/2  
          else:  
             return (1-s)/2  
291    
292    
293       def update(self,dt):       def update(self,dt):
# Line 128  class LevelSet(object): Line 304  class LevelSet(object):
304              self.__fcpde.setValue(Y = self.__phi-dt*inner(self.velocity,grad(phi_half)))              self.__fcpde.setValue(Y = self.__phi-dt*inner(self.velocity,grad(phi_half)))
305              self.__phi= self.__fcpde.getSolution()              self.__phi= self.__fcpde.getSolution()
306           self.__update_count += 1           self.__update_count += 1
307           if self.__update_count%self.__reinitialization_each == 0:           if self.__update_count%self.__reinit_each == 0:
308              self.__phi=self.__reinitialise(self.__phi)              self.__phi=self.__reinitialise(self.__phi)
309              if self.__FC:              if self.__FC:
310                  self.__fc.setInitialSolution(self.__phi+self.__diam)                  self.__fc.setInitialSolution(self.__phi+self.__diam)
# Line 153  class LevelSet(object): Line 329  class LevelSet(object):
329           # print "step size: dt (pos)= ",dtau           # print "step size: dt (pos)= ",dtau
330           # print "phi_p range:",inf(phi_p), sup(phi_p)           # print "phi_p range:",inf(phi_p), sup(phi_p)
331           # iter=0           # iter=0
332           # while (iter<=self.__reinitialization_steps_max):           # while (iter<=self.__reinit_max):
333           # phi_p=self.__reinitfc.solve(dtau)           # phi_p=self.__reinitfc.solve(dtau)
334           # print "phi_p range:",inf(phi_p), sup(phi_p)           # print "phi_p range:",inf(phi_p), sup(phi_p)
335           # iter+=1           # iter+=1
# Line 166  class LevelSet(object): Line 342  class LevelSet(object):
342           # print "step size: dt (neg)= ",dtau           # print "step size: dt (neg)= ",dtau
343           # print "phi_n range:",inf(phi_n), sup(phi_n)           # print "phi_n range:",inf(phi_n), sup(phi_n)
344           # iter=0           # iter=0
345           # while (iter<=self.__reinitialization_steps_max):           # while (iter<=self.__reinit_max):
346           # phi_n=self.__reinitfc.solve(dtau)           # phi_n=self.__reinitfc.solve(dtau)
347           # print "phi_n range:",inf(phi_n), sup(phi_n)           # print "phi_n range:",inf(phi_n), sup(phi_n)
348           # iter+=1           # iter+=1
# Line 187  class LevelSet(object): Line 363  class LevelSet(object):
363           iter =0           iter =0
364           # self.__reinitPde.setValue(q=whereNegative(abs(phi)-2*self.__h), r=phi)           # self.__reinitPde.setValue(q=whereNegative(abs(phi)-2*self.__h), r=phi)
365           # self.__reinitPde.setValue(r=phi)           # self.__reinitPde.setValue(r=phi)
366           while (iter<=self.__reinitialization_steps_max):           while (iter<=self.__reinit_max):
367                   phi_old=phi                   phi_old=phi
368                   if self.__reinitFC:                   if self.__reinitFC:
369                     phi = self.__reinitfc.solve(dtau)-self.__diam                     phi = self.__reinitfc.solve(dtau)-self.__diam
# Line 269  class LevelSet(object): Line 445  class LevelSet(object):
445          # saveVTK("test.%s.xml"%c,l=length(grad(self.__phi,fs))-1,s=s,phi=self.__phi)          # saveVTK("test.%s.xml"%c,l=length(grad(self.__phi,fs))-1,s=s,phi=self.__phi)
446    
447          dtau=f*inf(h/abs(s))          dtau=f*inf(h/abs(s))
448          while c < self.__reinitialization_steps_max: # and abs(diff) >= 0.01:          while c < self.__reinit_max: # and abs(diff) >= 0.01:
449            #            #
450            grad_phi=grad(self.__phi,fs)            grad_phi=grad(self.__phi,fs)
451            len_grad_phi=length(grad_phi)            len_grad_phi=length(grad_phi)

Legend:
Removed from v.1716  
changed lines
  Added in v.1788

  ViewVC Help
Powered by ViewVC 1.1.26