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

revision 1811 by ksteube, Thu Sep 25 23:11:13 2008 UTC revision 2158 by caltinay, Mon Dec 15 07:17:47 2008 UTC
# Line 36  import sys Line 36  import sys
36
37  class LevelSet:  class LevelSet:
38    """    """
39    The level set method tracking an interface defined by the zero contour of the level set function phi.    The level set method tracking an interface defined by the zero contour of the
40    it is assumed the phi(x)<0 defines the volume of interest.    level set function phi.
41      It is assumed that phi(x)<0 defines the volume of interest.
42
43    """    """
44    def __init__(self,domain,phi,reinit_max=10,reinit_each=2,smooth=2.):    def __init__(self,domain,phi,reinit_max=10,reinit_each=2,smooth=2.):
# Line 65  class LevelSet: Line 66  class LevelSet:
66      self.__h = inf(domain.getSize())      self.__h = inf(domain.getSize())
67      self.__smooth = smooth      self.__smooth = smooth
68      self.__n_step=0      self.__n_step=0
69
70    def __advect(self, velocity, dt):    def __advect(self, velocity, dt):
71      """      """
72      advects the level set function in the presense of a velocity field.      advects the level set function in the presence of a velocity field.
73
74      This implementation uses the 2-step Taylor-Galerkin method      This implementation uses the 2-step Taylor-Galerkin method
75      @param velocity: velocity field      @param velocity: velocity field
76      @param dt: dime increment      @param dt: time increment
77      @return: the advected level set function      @return: the advected level set function
78      """      """
79      Y = self.__phi-dt/2.*inner(velocity,grad(self.__phi))      Y = self.__phi-dt/2.*inner(velocity,grad(self.__phi))
80      self.__PDE.setValue(Y=Y)          self.__PDE.setValue(Y=Y)
81      phi_half = self.__PDE.getSolution()      phi_half = self.__PDE.getSolution()
82      Y = self.__phi-dt*inner(velocity,grad(phi_half))      Y = self.__phi-dt*inner(velocity,grad(phi_half))
83      self.__PDE.setValue(Y=Y)          self.__PDE.setValue(Y=Y)
84      phi = self.__PDE.getSolution()      phi = self.__PDE.getSolution()
85      print "LevelSet: Advection step done"      print "LevelSet: Advection step done"
86      return phi      return phi
87
88    def __reinitialise(self):    def __reinitialise(self):
89      """      """
90      reinializes the level set      reinitializes the level set
91
92      It solves the      It solves the...
93
94      @return: reinitalized level set      @return: reinitialized level set
95      """      """
96      phi=self.__phi      phi=self.__phi
97      s = sign(phi.interpolate(Function(self.__domain)))      s = sign(phi.interpolate(Function(self.__domain)))
# Line 116  class LevelSet: Line 117  class LevelSet:
117
118    def getTimeStepSize(self,velocity):    def getTimeStepSize(self,velocity):
119         """         """
120         returns a new dt for a given velocity using the courant coundition         returns a new dt for a given velocity using the courant condition
121
122         @param velocity: velocity field         @param velocity: velocity field
123         """         """
124         self.__velocity=velocity         self.__velocity=velocity
125         dt=0.5*self.__h/sup(length(velocity))         dt=0.5*self.__h/sup(length(velocity))
126         return dt         return dt
127
128    def update(self,dt):    def update(self,dt):
129        """        """
130        sets a new velocity and updates the level set fuction        sets a new velocity and updates the level set function
131
132        @param dt: time step forward        @param dt: time step forward
133        """        """
# Line 137  class LevelSet: Line 138  class LevelSet:
138
139    def update_phi(self, velocity, dt):    def update_phi(self, velocity, dt):
140        """        """
141        updates phi under the presense of a velocity field        updates phi under the presence of a velocity field
142
143          If dt is small this call is equivalent to call
144
If dt is small this call is equivalent to call

145        dt=LevelSet.getTimeStepSize(velocity)        dt=LevelSet.getTimeStepSize(velocity)
146        phi=LevelSet.update(dt)        phi=LevelSet.update(dt)
147
148        otherwise substepping is used.        otherwise substepping is used.
149        @param velocity: velocity field        @param velocity: velocity field
150        @param dt: time step forward        @param dt: time step forward
# Line 159  class LevelSet: Line 160  class LevelSet:
160
161    def getVolume(self):    def getVolume(self):
162      """      """
163      return the volume of the phi(x)<0 region      returns the volume of the phi(x)<0 region
164      """      """
165      return integrate(whereNegative(self.__phi.interpolate(Function(self.__domain))))      return integrate(whereNegative(self.__phi.interpolate(Function(self.__domain))))
166
167    def getSurface(self,rel_width_factor=0.5):    def getSurface(self,rel_width_factor=0.5):
168      """      """
169      return a mask for phi(x)=1 region      returns a mask for the phi(x)=1 region
170
171      @param rel_width_factor: relative wideth of region around zero contour.      @param rel_width_factor: relative width of region around zero contour.
172      """      """
173      return whereNegative(abs(self.__phi)-rel_width_factor*self.__h)      return whereNegative(abs(self.__phi)-rel_width_factor*self.__h)
174
175    def getH(self):    def getH(self):
176       """       """
177       returns mesh size       returns the mesh size
178       """       """
179       return self.__h       return self.__h
180
181    def getDomain(self):    def getDomain(self):
182       """       """
183       returns domain       returns the domain
184       """       """
185       return self.__domain       return self.__domain
186
# Line 191  class LevelSet: Line 192  class LevelSet:
192
193    def update_parameter_sharp(self, param_neg=-1, param_pos=1, phi=None):    def update_parameter_sharp(self, param_neg=-1, param_pos=1, phi=None):
194      """      """
195      creates a function whith param_neg where phi<0 and param_pos where phi>0 (no smoothing)      creates a function with param_neg where phi<0 and param_pos where phi>0
196        (no smoothing)
197
198      @param param_neg: value of parameter on the negative side (phi<0)      @param param_neg: value of parameter on the negative side (phi<0)
199      @param param_pos: value of parameter on the positve side (phi>0)      @param param_pos: value of parameter on the positive side (phi>0)
200      @param phi: level set funtion to be used. if not present the current level set is used.      @param phi: level set function to be used. if not present the current level set is used.
201      """      """
202      mask_neg = whereNegative(self.__phi)      mask_neg = whereNegative(self.__phi)
203      mask_pos = whereNonNegative(self.__phi)      mask_pos = whereNonNegative(self.__phi)
# Line 204  class LevelSet: Line 206  class LevelSet:
206
207    def update_parameter(self, param_neg=-1, param_pos=1, phi=None, smoothing_width=None):    def update_parameter(self, param_neg=-1, param_pos=1, phi=None, smoothing_width=None):
208      """      """
209      creates a smoothed function whith param_neg where phi<0 and param_pos where phi>0 which is smoothed over a length      creates a smoothed function with param_neg where phi<0 and param_pos where
210      smoothing_width accross the interface      phi>0 which is smoothed over a length smoothing_width across the interface
211
212      @param smoothing_width: width of the smoothing zone relative to mesh size. If not present the initial value of C{smooth} is used.      @param smoothing_width: width of the smoothing zone relative to mesh size.
213                                If not present the initial value of C{smooth} is used.
214      """      """
215      if smoothing_width==None: smoothing_width = self.__smooth      if smoothing_width==None: smoothing_width = self.__smooth
216      if phi==None: phi=self.__phi      if phi==None: phi=self.__phi
# Line 216  class LevelSet: Line 219  class LevelSet:
219
220    def __makeInterface(self,phi,smoothing_width):    def __makeInterface(self,phi,smoothing_width):
221        """        """
222        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        creates a smooth interface from -1 to 1 over the length
223        and 1 where the level set is 1        2*h*smoothing_width where -1 is used where the level set is negative and
224          1 where the level set is 1
225        """        """
226        s=smoothing_width*self.__h        s=smoothing_width*self.__h
227        phi_on_h=interpolate(phi,Function(self.__domain))                phi_on_h=interpolate(phi,Function(self.__domain))
228        mask_neg = whereNonNegative(-s-phi_on_h)        mask_neg = whereNonNegative(-s-phi_on_h)
229        mask_pos = whereNonNegative(phi_on_h-s)        mask_pos = whereNonNegative(phi_on_h-s)
# Line 229  class LevelSet: Line 233  class LevelSet:
233
234    def makeCharacteristicFunction(self, contour=0, phi=None, positiveSide=True, smoothing_width=None):    def makeCharacteristicFunction(self, contour=0, phi=None, positiveSide=True, smoothing_width=None):
235        """        """
236        makes a smooth charateristic function of the region phi(x)>contour if positiveSide and phi(x)<contour otherwise.        makes a smooth characteristic function of the region phi(x)>contour if
237          positiveSide and phi(x)<contour otherwise.
238
239        @param phi: level set funtion to be used. if not present the current level set is used.        @param phi: level set function to be used. If not present the current level set is used.
240        @param smoothing_width: width of the smoothing zone relative to mesh size. If not present the initial value of C{smooth} is used.        @param smoothing_width: width of the smoothing zone relative to mesh size.
241                                  If not present the initial value of C{smooth} is used.
242        """        """
243        if phi==None: phi=self.__phi        if phi==None: phi=self.__phi
244        if smoothing_width == None: smoothing_width=self.__smooth        if smoothing_width == None: smoothing_width=self.__smooth
# Line 252  class LevelSet2(object): Line 258  class LevelSet2(object):
258           """           """
259           initialize model           initialize model
260           """           """
261           self.__domain = phi.getDomain()           self.__domain = phi.getDomain()
262           x=self.__domain.getX()           x=self.__domain.getX()
263           diam=0           diam=0
264           for i in range(self.__domain.getDim()):           for i in range(self.__domain.getDim()):
# Line 276  class LevelSet2(object): Line 282  class LevelSet2(object):
282              self.__fc.setInitialSolution(phi+self.__diam)              self.__fc.setInitialSolution(phi+self.__diam)
283           else:           else:
284              self.__fcpde = LinearPDE(self.__domain,numEquations=1, numSolutions=1)              self.__fcpde = LinearPDE(self.__domain,numEquations=1, numSolutions=1)
285              self.__fcpde.setSolverMethod(solver=LinearPDE.LUMPING)                    self.__fcpde.setSolverMethod(solver=LinearPDE.LUMPING)
286              self.__fcpde.setReducedOrderOn()              self.__fcpde.setReducedOrderOn()
287              self.__fcpde.setSymmetryOn()              self.__fcpde.setSymmetryOn()
288              self.__fcpde.setValue(D=1.)              self.__fcpde.setValue(D=1.)
# Line 289  class LevelSet2(object): Line 295  class LevelSet2(object):
295              self.__reinitfc.setTolerance(1.e-5)              self.__reinitfc.setTolerance(1.e-5)
296           else:           else:
297              self.__reinitPde = LinearPDE(self.__domain,numEquations=1, numSolutions=1)              self.__reinitPde = LinearPDE(self.__domain,numEquations=1, numSolutions=1)
298              # self.__reinitPde.setSolverMethod(solver=LinearPDE.LUMPING)                    # self.__reinitPde.setSolverMethod(solver=LinearPDE.LUMPING)
299              self.__reinitPde.setReducedOrderOn()              self.__reinitPde.setReducedOrderOn()
300              self.__reinitPde.setSymmetryOn()              self.__reinitPde.setSymmetryOn()
301              self.__reinitPde.setValue(D=1.)              self.__reinitPde.setValue(D=1.)
# Line 309  class LevelSet2(object): Line 315  class LevelSet2(object):
315           self.__smoothed_char=self.__makeInterface(self.__phi)           self.__smoothed_char=self.__makeInterface(self.__phi)
316
317

318       def update(self,dt):       def update(self,dt):
319           """           """
320           sets a new velocity and updates the level set fuction           sets a new velocity and updates the level set function
321
322           @param dt: time step forward           @param dt: time step forward
323           """           """
# Line 358  class LevelSet2(object): Line 362  class LevelSet2(object):
362           # phi_n=-whereNegative(phi)*phi           # phi_n=-whereNegative(phi)*phi
363           # self.__reinitfc.setInitialSolution(phi_n)           # self.__reinitfc.setInitialSolution(phi_n)
364           # self.__reinitfc.setValue(C=-whereNegative(s)*w,Y=-whereNegative(s)*s,q=whereNonNegative(phi))           # self.__reinitfc.setValue(C=-whereNegative(s)*w,Y=-whereNegative(s)*s,q=whereNonNegative(phi))
365           # # dtau=self.__reinitfc.getSafeTimeStepSize()           # # dtau=self.__reinitfc.getSafeTimeStepSize()
366           # dtau=self.__h           # dtau=self.__h
367           # print "step size: dt (neg)= ",dtau           # print "step size: dt (neg)= ",dtau
368           # print "phi_n range:",inf(phi_n), sup(phi_n)           # print "phi_n range:",inf(phi_n), sup(phi_n)
# Line 376  class LevelSet2(object): Line 380  class LevelSet2(object):
380           if self.__reinitFC:           if self.__reinitFC:
381               self.__reinitfc.setValue(C=-w,Y=s)               self.__reinitfc.setValue(C=-w,Y=s)
382               self.__reinitfc.setInitialSolution(phi+self.__diam)               self.__reinitfc.setInitialSolution(phi+self.__diam)
383               dtau=self.__reinitfc.getSafeTimeStepSize()               dtau=self.__reinitfc.getSafeTimeStepSize()
384               # dtau = 0.5*inf(Function(self.__domain).getSize())               # dtau = 0.5*inf(Function(self.__domain).getSize())
385           else:           else:
386               dtau = 0.5*inf(Function(self.__domain).getSize())               dtau = 0.5*inf(Function(self.__domain).getSize())
# Line 408  class LevelSet2(object): Line 412  class LevelSet2(object):
412           return out           return out
413
414
415       #================ things from here onwards are not used nor tested: ==========================================       #=========== things from here onwards are not used nor tested: ===========
416
417
418       def getCharacteristicFunction(self):       def getCharacteristicFunction(self):
419           return self.__smoothed_char           return self.__smoothed_char
420
421

422       def RK2(self,L):       def RK2(self,L):
423             k0=L(phi)             k0=L(phi)
424             phi_1=phi+dt*k0             phi_1=phi+dt*k0
425             k1=L(phi_1)             k1=L(phi_1)
426             phi=phi+dt/2*(k0+k1)             phi=phi+dt/2*(k0+k1)
427
428       def RK3(self,L):       def RK3(self,L):
429             k0=L(phi)             k0=L(phi)
430             phi_1=phi+dt*L(phi)             phi_1=phi+dt*L(phi)
# Line 429  class LevelSet2(object): Line 432  class LevelSet2(object):
432             phi_2=phi+dt/4*(k0+k1)             phi_2=phi+dt/4*(k0+k1)
433             k2=L(phi_2)             k2=L(phi_2)
434             phi=phi+dt/6*(k0+4*k2+K1)             phi=phi+dt/6*(k0+4*k2+K1)
435
436       def TG(self,L):       def TG(self,L):
437             k0=L(phi)             k0=L(phi)
438             phi_1=phi+dt/2*k0             phi_1=phi+dt/2*k0
439             k1=L(phi_1)             k1=L(phi_1)
440             phi=phi+dt*k1             phi=phi+dt*k1
441
442
443       def __reinitialise_old(self):       def __reinitialise_old(self):
444          #=============================================          #=============================================
445          f=0.1          f=0.1
# Line 474  class LevelSet2(object): Line 478  class LevelSet2(object):
478            # s=self.__makeInterface(1.)            # s=self.__makeInterface(1.)
480            # s=s*s2            # s=s*s2
481            # phi_on_h=interpolate(self.__phi,fs)                    # phi_on_h=interpolate(self.__phi,fs)
482            # self.__reinitPde.setValue(Y =self.__phi+dtau/2*s*(1.-len_grad_phi))            # self.__reinitPde.setValue(Y =self.__phi+dtau/2*s*(1.-len_grad_phi))
483            # phi_half=self.__reinitPde.getSolution()            # phi_half=self.__reinitPde.getSolution()
484            # self.__reinitPde.setValue(Y =self.__phi+dtau*s*(1.-length(grad(phi_half,fs))))            # self.__reinitPde.setValue(Y =self.__phi+dtau*s*(1.-length(grad(phi_half,fs))))
# Line 503  class LevelSet2(object): Line 507  class LevelSet2(object):
507          s=abs(self. __makeInterface(2.))          s=abs(self. __makeInterface(2.))
510
511          #----          #----
512          # aphi_on_h=abs(interpolate(self.__phi,self.__h.getFunctionSpace()))          # aphi_on_h=abs(interpolate(self.__phi,self.__h.getFunctionSpace()))
513          # q=2*self.__h          # q=2*self.__h
# Line 535  class LevelSet2(object): Line 539  class LevelSet2(object):
539
540       def getVolumeOfNegativeDomain(self):       def getVolumeOfNegativeDomain(self):
541           """           """
542           return the current volume of domain with phi<0.           returns the current volume of domain with phi<0.
543           """           """
544           return integrate((1.-self.__makeInterface(1.))/2.)           return integrate((1.-self.__makeInterface(1.))/2.)
545
546
547       def getBoundingBoxOfNegativeDomain(self):       def getBoundingBoxOfNegativeDomain(self):
548           """           """
549           get the height of the region with  phi<0           returns the height of the region with phi<0
550           """           """
551           fs=self.__h.getFunctionSpace()           fs=self.__h.getFunctionSpace()