/[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 3070 by jfenwick, Thu Jan 28 02:03:15 2010 UTC revision 3071 by gross, Wed Jul 21 05:37:30 2010 UTC
# Line 20  http://www.opensource.org/licenses/osl-3 Line 20  http://www.opensource.org/licenses/osl-3
20  __url__="https://launchpad.net/escript-finley"  __url__="https://launchpad.net/escript-finley"
21    
22  from esys.escript import *  from esys.escript import *
23  from esys.escript.linearPDEs import LinearPDE, TransportPDE  from esys.escript.linearPDEs import LinearPDE, SingleTransportPDE
24  from esys.escript.pdetools import Projector  from esys.escript.pdetools import Projector
25  from esys import finley  from esys import finley
26  import math  import math
27    
28    USE_OLD_VERSION=False
29    USE_OLD_VERSION_REINIT=False
 from esys.escript import *  
 from esys.finley import finley  
 from esys.escript.linearPDEs import LinearPDE  
 from esys.escript.pdetools import Projector  
 import sys  
30    
31    
32  class LevelSet:  class LevelSet:
33    """    """
34    The level set method tracking an interface defined by the zero contour of the    The level set method tracking an interface defined by the zero contour of the
35    level set function phi.    level set function phi which defines the signed distance of a point x from the
36      interface. The contour phi(x)=0 defines the interface.
   It is assumed that phi(x)<0 defines the volume of interest.  
37    
38      It is assumed that phi(x)<0 defines the volume of interest,
39      phi(x)>0 the outside world.
40    """    """
41    def __init__(self,domain,phi,reinit_max=10,reinit_each=2,smooth=2.):    def __init__(self,phi,reinit_max=10,reinitialize_after=20,smooth=2., useReducedOrder=False):
42      """      """
43      Sets up the level set method.      Sets up the level set method.
44    
45      :param domain: the domain where the level set is used      :param domain: the domain where the level set is used
46      :param phi: the initial level set function      :param phi: the initial level set function
47      :param reinit_max: maximum number of reinitialization steps      :param reinit_max: maximum number of reinitialization steps
48      :param reinit_each: ``phi`` is reinitialized every ``reinit_each`` step      :param reinit_after: ``phi`` is reinitialized every ``reinit_after`` step
49      :param smooth: smoothing width      :param smooth: smoothing width
50      """      """
51      self.__domain = domain      self.__domain = phi.getDomain()
52      self.__phi = phi      self.__phi = phi
53    
54        if USE_OLD_VERSION:
55           self.__PDE = LinearPDE(self.__domain)
56           if useReducedOrder: self.__PDE.setReducedOrderOn()
57           self.__PDE.setSolverMethod(solver=LinearPDE.PCG)
58           self.__PDE.setValue(D=1.0)
59        else:
60           self.__PDE=SingleTransportPDE(self.__domain)
61           if useReducedOrder: self.__PDE.setReducedOrderOn()
62           self.__PDE.setValue(M=1.0)
63    
64        if USE_OLD_VERSION_REINIT:
65           self.__reinitPDE = LinearPDE(self.__domain, numEquations=1)
66           self.__reinitPDE.getSolverOption().setSolverMethod(solver=LinearPDE.LUMPING)
67           if useReducedOrder: self.__reinitPDE.setReducedOrderOn()
68           self.__reinitPDE.setValue(D=1.0)
69        else:
70           self.__reinitPDE=SingleTransportPDE(self.__domain)
71           if useReducedOrder: self.__reinitPDE.setReducedOrderOn()
72           self.__reinitPDE.setValue(M=1.0)
73    
74        # revise:
75      self.__reinit_max = reinit_max      self.__reinit_max = reinit_max
76      self.__reinit_each = reinit_each      self.__reinit_after = reinitialize_after
77      self.__PDE = LinearPDE(domain)      self.__h = inf(self.__domain.getSize())
     self.__PDE.setReducedOrderOn()  
     self.__PDE.setValue(D=1.0)  
     self.__PDE.setSolverMethod(solver=LinearPDE.PCG)  
     self.__reinitPDE = LinearPDE(domain, numEquations=1)  
     self.__reinitPDE.setReducedOrderOn()  
     self.__reinitPDE.setValue(D=1.0)  
     self.__reinitPDE.setSolverMethod(solver=LinearPDE.LUMPING)  
     self.__h = inf(domain.getSize())  
78      self.__smooth = smooth      self.__smooth = smooth
79      self.__n_step=0      self.__n_step=0
80    
81    def __advect(self, velocity, dt):    def getH(self):
82      """        """
83      Advects the level set function in the presence of a velocity field.        Returns the mesh size.
84          """
85      This implementation uses the 2-step Taylor-Galerkin method.        return self.__h
86    
87      :param velocity: velocity field    def getDomain(self):
88      :param dt: time increment        """
89      :return: the advected level set function        Returns the domain.
90      """        """
91      Y = self.__phi-dt/2.*inner(velocity,grad(self.__phi))        return self.__domain
     self.__PDE.setValue(Y=Y)  
     phi_half = self.__PDE.getSolution()  
     Y = self.__phi-dt*inner(velocity,grad(phi_half))  
     self.__PDE.setValue(Y=Y)  
     phi = self.__PDE.getSolution()  
     print "LevelSet: Advection step done"  
     return phi  
92    
93    def __reinitialise(self):    def getAdvectionSolverOptions(self):
94      """        """
95      Reinitializes the level set.        Returns the solver options for the interface advective.
96          """
97          return self.__PDE.getSolverOptions()
98    
99      It solves the PDE...    def getReinitializationSolverOptions(self):
100          """
101          Returns the options of the solver for the reinitialization
102          """
103          return self.__reinitPDE.getSolverOption()
104    
105      :return: reinitialized level set    def getLevelSetFunction(self):
106      """        """
107      phi=self.__phi        Returns the level set function.
108      s = sign(phi.interpolate(Function(self.__domain)))        """
109      g=grad(phi)        return self.__phi
     w = s*g/(length(g)+1e-10)  
     dtau = 0.5*self.__h  
     iter =0  
     mask = whereNegative(abs(phi)-1.*self.__h)  
     self.__reinitPDE.setValue(q=mask, r=phi)  
     g=grad(phi)  
     while (iter<=self.__reinit_max):  
       Y = phi+(dtau/2.)*(s-inner(w,g))  
       self.__reinitPDE.setValue(Y = Y)  
       phi_half = self.__reinitPDE.getSolution()  
       Y = phi+dtau*(s-inner(w,grad(phi_half)))  
       self.__reinitPDE.setValue(Y = Y)  
       phi = self.__reinitPDE.getSolution()  
       g=grad(phi)  
       error = Lsup(length(g)*whereNegative(abs(phi.interpolate(g.getFunctionSpace()))-3.0*self.__h))  
       print "LevelSet:reinitialization: iteration :", iter, " error:", error  
       iter +=1  
     return phi  
110    
111    def getTimeStepSize(self,velocity):    def getTimeStepSize(self,velocity):
112         """         """
# Line 123  class LevelSet: Line 114  class LevelSet:
114    
115         :param velocity: velocity field         :param velocity: velocity field
116         """         """
117         self.__velocity=velocity         if USE_OLD_VERSION:
118         dt=0.5*self.__h/sup(length(velocity))            self.__velocity=velocity
119              dt=0.5*self.__h/sup(length(velocity))
120           else:
121              self.__PDE.setValue(C=-velocity)
122              dt=self.__PDE.getSafeTimeStepSize()
123    
124         return dt         return dt
125    
126    def update(self,dt):    def update(self,dt):
# Line 133  class LevelSet: Line 129  class LevelSet:
129    
130        :param dt: time step forward        :param dt: time step forward
131        """        """
132        self.__phi=self.__advect(self.__velocity, dt)        self.__phi=self.__advect(dt)
133        self.__n_step+=1        self.__n_step+=1
134        if self.__n_step%self.__reinit_each ==0: self.__phi = self.__reinitialise()        if self.__n_step%self.__reinit_after ==0: self.__phi = self.__reinitialise()
135        return self.__phi        return self.__phi
136    
137    def update_phi(self, velocity, dt):    def update_phi(self, velocity, dt):
# Line 153  class LevelSet: Line 149  class LevelSet:
149        :param dt: time step forward        :param dt: time step forward
150        """        """
151        dt2=self.getTimeStepSize(velocity)        dt2=self.getTimeStepSize(velocity)
152        n=math.ceil(dt/dt2)        n=int(math.ceil(dt/dt2))
153        dt_new=dt/n        dt_new=dt/n
154        for i in range(n):        for i in range(n):
155             phi=self.update(dt_new)             phi=self.update(dt_new)
            t+=dt_new  
156        return phi        return phi
157    
158      def __advect(self,  dt):
159        """
160        Advects the level set function in the presence of a velocity field.
161    
162    def getVolume(self):      :param dt: time increment
163        """      :return: the advected level set function
164        Returns the volume of the *phi(x)<0* region.      """
165        """      if USE_OLD_VERSION:
166        return integrate(whereNegative(self.__phi.interpolate(Function(self.__domain))))         velocity=self.__velocity
167           Y = self.__phi-dt/2.*inner(velocity,grad(self.__phi))
168           self.__PDE.setValue(Y=Y)
169           phi_half = self.__PDE.getSolution()
170           Y = self.__phi-dt*inner(velocity,grad(phi_half))
171           self.__PDE.setValue(Y=Y)
172           phi = self.__PDE.getSolution()
173        else:
174           phi = self.__PDE.getSolution(dt, self.__phi)
175        print "LevelSet: Advection step done"
176        return phi
177    
178    def getSurface(self,rel_width_factor=0.5):    #==============================================================================================
179        """    def __reinitialise(self):
180        Returns a mask for the *phi(x)=1* region      """
181        Reinitializes the level set.
182    
183        :param rel_width_factor: relative width of region around zero contour      It solves the PDE...
       """  
       return whereNegative(abs(self.__phi)-rel_width_factor*self.__h)  
184    
185    def getH(self):      :return: reinitialized level set
186        """      """
187        Returns the mesh size.      phi=self.__phi
188        """      s = sign(phi.interpolate(Function(self.__domain)))
189        return self.__h      g=grad(phi)
190        w = s*g/(length(g)+1e-10)
191        dtau = 0.5*self.__h
192        iter =0
193        mask = whereNegative(abs(phi)-1.*self.__h)
194        self.__reinitPDE.setValue(q=mask, r=phi)
195        g=grad(phi)
196        while (iter<=self.__reinit_max):
197          Y = phi+(dtau/2.)*(s-inner(w,g))
198          self.__reinitPDE.setValue(Y = Y)
199          phi_half = self.__reinitPDE.getSolution()
200          Y = phi+dtau*(s-inner(w,grad(phi_half)))
201          self.__reinitPDE.setValue(Y = Y)
202          phi = self.__reinitPDE.getSolution()
203          g=grad(phi)
204          error = Lsup(length(g)*whereNegative(abs(phi.interpolate(g.getFunctionSpace()))-3.0*self.__h))
205          print "LevelSet:reinitialization: iteration :", iter, " error:", error
206          iter +=1
207        return phi
208    
   def getDomain(self):  
       """  
       Returns the domain.  
       """  
       return self.__domain  
209    
210    def getLevelSetFunction(self):    def getVolume(self):
211        """        """
212        Returns the level set function.        Returns the volume of the *phi(x)<0* region.
213        """        """
214        return self.__phi        return integrate(whereNegative(self.__phi.interpolate(Function(self.__domain))))
215    
216    
217    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):
218        """        """
# Line 256  class LevelSet: Line 277  class LevelSet:
277        else:        else:
278            return (1-s)/2            return (1-s)/2
279    
   def setTolerance(self,tolerance=1e-3):  
     self.__PDE.setTolerance(tolerance)  
     self.__reinitPDE.setTolerance(tolerance)  
   
280    
281  class LevelSet2(object):  class LevelSet2(object):
282       def __init__(self,phi,reinit_max=10,reinit_each=2,smooth=2.):       def __init__(self,phi,reinit_max=10,reinit_after=2,smooth=2.):
283           """           """
284           Initializes the model.           Initializes the model.
285           """           """
# Line 274  class LevelSet2(object): Line 291  class LevelSet2(object):
291               diam+=(inf(xi)-sup(xi))**2               diam+=(inf(xi)-sup(xi))**2
292           self.__diam=sqrt(diam)           self.__diam=sqrt(diam)
293           self.__h = sup(Function(self.__domain).getSize())           self.__h = sup(Function(self.__domain).getSize())
294           self.__reinit_each=reinit_each           self.__reinit_after=reinit_after
295           self.__reinit_max = reinit_max           self.__reinit_max = reinit_max
296           self.__smooth = smooth           self.__smooth = smooth
297           self.__phi = phi           self.__phi = phi
# Line 337  class LevelSet2(object): Line 354  class LevelSet2(object):
354              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)))
355              self.__phi= self.__fcpde.getSolution()              self.__phi= self.__fcpde.getSolution()
356           self.__update_count += 1           self.__update_count += 1
357           if self.__update_count%self.__reinit_each == 0:           if self.__update_count%self.__reinit_after == 0:
358              self.__phi=self.__reinitialise(self.__phi)              self.__phi=self.__reinitialise(self.__phi)
359              if self.__FC:              if self.__FC:
360                  self.__fc.setInitialSolution(self.__phi+self.__diam)                  self.__fc.setInitialSolution(self.__phi+self.__diam)

Legend:
Removed from v.3070  
changed lines
  Added in v.3071

  ViewVC Help
Powered by ViewVC 1.1.26