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

revision 2158 by caltinay, Mon Dec 15 07:17:47 2008 UTC revision 2169 by caltinay, Wed Dec 17 03:08:58 2008 UTC
# Line 38  class LevelSet: Line 38  class LevelSet:
38    """    """
39    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
40    level set function phi.    level set function phi.
41
42    It is assumed that phi(x)<0 defines the volume of interest.    It is assumed that phi(x)<0 defines the volume of interest.
43
44    """    """
45    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.):
46      """      """
47      set up the level set method      Sets up the level set method.
48
49      @param domain: the domain where the level set is used      @param domain: the domain where the level set is used
50      @param phi: the initial level set function      @param phi: the initial level set function
51      @param reinit_max: maximum number of reinitalization steps      @param reinit_max: maximum number of reinitialization steps
52      @param reinit_each: phi is reinitialized every reinit_each step      @param reinit_each: C{phi} is reinitialized every C{reinit_each} step
53      @param smooth: smoothing width      @param smooth: smoothing width
54      """      """
55      self.__domain = domain      self.__domain = domain
# Line 69  class LevelSet: Line 70  class LevelSet:
70
72      """      """
73      advects the level set function in the presence of a velocity field.      Advects the level set function in the presence of a velocity field.
74
75        This implementation uses the 2-step Taylor-Galerkin method.
76
This implementation uses the 2-step Taylor-Galerkin method
77      @param velocity: velocity field      @param velocity: velocity field
78      @param dt: time increment      @param dt: time increment
79      @return: the advected level set function      @return: the advected level set function
# Line 87  class LevelSet: Line 89  class LevelSet:
89
90    def __reinitialise(self):    def __reinitialise(self):
91      """      """
92      reinitializes the level set      Reinitializes the level set.
93
94      It solves the...      It solves the PDE...
95
96      @return: reinitialized level set      @return: reinitialized level set
97      """      """
# Line 117  class LevelSet: Line 119  class LevelSet:
119
120    def getTimeStepSize(self,velocity):    def getTimeStepSize(self,velocity):
121         """         """
122         returns a new dt for a given velocity using the courant condition         Returns a new C{dt} for a given C{velocity} using the Courant condition.
123
124         @param velocity: velocity field         @param velocity: velocity field
125         """         """
# Line 127  class LevelSet: Line 129  class LevelSet:
129
130    def update(self,dt):    def update(self,dt):
131        """        """
132        sets a new velocity and updates the level set function        Sets a new velocity and updates the level set function.
133
134        @param dt: time step forward        @param dt: time step forward
135        """        """
# Line 138  class LevelSet: Line 140  class LevelSet:
140
141    def update_phi(self, velocity, dt):    def update_phi(self, velocity, dt):
142        """        """
143        updates phi under the presence of a velocity field        Updates C{phi} under the presence of a velocity field.
144
145        If dt is small this call is equivalent to call        If dt is small this call is equivalent to::
146
147        dt=LevelSet.getTimeStepSize(velocity)            dt=LevelSet.getTimeStepSize(velocity)
148        phi=LevelSet.update(dt)            phi=LevelSet.update(dt)
149
150        otherwise substepping is used.        otherwise substepping is used.
151
152        @param velocity: velocity field        @param velocity: velocity field
153        @param dt: time step forward        @param dt: time step forward
154        """        """
# Line 159  class LevelSet: Line 162  class LevelSet:
162
163
164    def getVolume(self):    def getVolume(self):
165      """        """
166      returns the volume of the phi(x)<0 region        Returns the volume of the M{phi(x)<0} region.
167      """        """
168      return integrate(whereNegative(self.__phi.interpolate(Function(self.__domain))))        return integrate(whereNegative(self.__phi.interpolate(Function(self.__domain))))
169
170    def getSurface(self,rel_width_factor=0.5):    def getSurface(self,rel_width_factor=0.5):
171      """        """
172      returns a mask for the phi(x)=1 region        Returns a mask for the M{phi(x)=1} region
173
174      @param rel_width_factor: relative width of region around zero contour.        @param rel_width_factor: relative width of region around zero contour
175      """        """
176      return whereNegative(abs(self.__phi)-rel_width_factor*self.__h)        return whereNegative(abs(self.__phi)-rel_width_factor*self.__h)
177
178    def getH(self):    def getH(self):
179       """        """
180       returns the mesh size        Returns the mesh size.
181       """        """
182       return self.__h        return self.__h
183
184    def getDomain(self):    def getDomain(self):
185       """        """
186       returns the domain        Returns the domain.
187       """        """
188       return self.__domain        return self.__domain
189
190    def getLevelSetFunction(self):    def getLevelSetFunction(self):
191        """        """
192        returns the level set function        Returns the level set function.
193        """        """
194        return self.__phi        return self.__phi
195
196    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):
197      """        """
198      creates a function with param_neg where phi<0 and param_pos where phi>0        Creates a function with C{param_neg} where C{phi<0} and C{param_pos}
199      (no smoothing)        where C{phi>0} (no smoothing).
200
201      @param param_neg: value of parameter on the negative side (phi<0)        @param param_neg: value of parameter on the negative side (phi<0)
202      @param param_pos: value of parameter on the positive side (phi>0)        @param param_pos: value of parameter on the positive side (phi>0)
203      @param phi: level set function 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
204      """                    level set is used.
209          return param
210
211    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):
212      """        """
213      creates a smoothed function with param_neg where phi<0 and param_pos where        Creates a smoothed function with C{param_neg} where C{phi<0} and
214      phi>0 which is smoothed over a length smoothing_width across the interface        C{param_pos} where C{phi>0} which is smoothed over a length
215          C{smoothing_width} across the interface.
216
217      @param smoothing_width: width of the smoothing zone relative to mesh size.        @param smoothing_width: width of the smoothing zone relative to mesh size.
218                              If not present the initial value of C{smooth} is used.                                If not present the initial value of C{smooth} is
219      """                                used.
220      if smoothing_width==None: smoothing_width = self.__smooth        """
221      if phi==None: phi=self.__phi        if smoothing_width==None: smoothing_width = self.__smooth
222      s=self.__makeInterface(phi,smoothing_width)        if phi==None: phi=self.__phi
223      return ((param_pos-param_neg)*s+param_pos+param_neg)/2        s=self.__makeInterface(phi,smoothing_width)
224          return ((param_pos-param_neg)*s+param_pos+param_neg)/2
225
226    def __makeInterface(self,phi,smoothing_width):    def __makeInterface(self,phi,smoothing_width):
227        """        """
228        creates a smooth interface from -1 to 1 over the length        Creates a smooth interface from -1 to 1 over the length
229        2*h*smoothing_width where -1 is used where the level set is negative and        M{2*h*smoothing_width} where -1 is used where the level set is negative
230        1 where the level set is 1        and 1 where the level set is 1.
231        """        """
232        s=smoothing_width*self.__h        s=smoothing_width*self.__h
233        phi_on_h=interpolate(phi,Function(self.__domain))        phi_on_h=interpolate(phi,Function(self.__domain))
# Line 233  class LevelSet: Line 239  class LevelSet:
239
240    def makeCharacteristicFunction(self, contour=0, phi=None, positiveSide=True, smoothing_width=None):    def makeCharacteristicFunction(self, contour=0, phi=None, positiveSide=True, smoothing_width=None):
241        """        """
242        makes a smooth characteristic function of the region phi(x)>contour if        Makes a smooth characteristic function of the region C{phi(x)>contour} if
243        positiveSide and phi(x)<contour otherwise.        C{positiveSide} and C{phi(x)<contour} otherwise.
244
245        @param phi: level set function 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
246                      level set is used.
247        @param smoothing_width: width of the smoothing zone relative to mesh size.        @param smoothing_width: width of the smoothing zone relative to mesh size.
248                                If not present the initial value of C{smooth} is used.                                If not present the initial value of C{smooth} is
249                                  used.
250        """        """
251        if phi==None: phi=self.__phi        if phi==None: phi=self.__phi
252        if smoothing_width == None: smoothing_width=self.__smooth        if smoothing_width == None: smoothing_width=self.__smooth
# Line 256  class LevelSet: Line 264  class LevelSet:
264  class LevelSet2(object):  class LevelSet2(object):
265       def __init__(self,phi,reinit_max=10,reinit_each=2,smooth=2.):       def __init__(self,phi,reinit_max=10,reinit_each=2,smooth=2.):
266           """           """
267           initialize model           Initializes the model.
268           """           """
269           self.__domain = phi.getDomain()           self.__domain = phi.getDomain()
270           x=self.__domain.getX()           x=self.__domain.getX()
# Line 317  class LevelSet2(object): Line 325  class LevelSet2(object):
325
326       def update(self,dt):       def update(self,dt):
327           """           """
328           sets a new velocity and updates the level set function           Sets a new velocity and updates the level set function.
329
330           @param dt: time step forward           @param dt: time step forward
331           """           """
# Line 539  class LevelSet2(object): Line 547  class LevelSet2(object):
547
548       def getVolumeOfNegativeDomain(self):       def getVolumeOfNegativeDomain(self):
549           """           """
550           returns the current volume of domain with phi<0.           Returns the current volume of domain with phi<0.
551           """           """
552           return integrate((1.-self.__makeInterface(1.))/2.)           return integrate((1.-self.__makeInterface(1.))/2.)
553
554
555       def getBoundingBoxOfNegativeDomain(self):       def getBoundingBoxOfNegativeDomain(self):
556           """           """
557           returns the height of the region with phi<0           Returns the height of the region with phi<0.
558           """           """
559           fs=self.__h.getFunctionSpace()           fs=self.__h.getFunctionSpace()