/[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 1473 by gross, Fri Apr 4 07:49:18 2008 UTC revision 1788 by gross, Mon Sep 15 02:20:20 2008 UTC
# Line 1  Line 1 
1  from esys.escript import *  from esys.escript import *
2  from esys.escript.linearPDEs import LinearPDE  from esys.escript.linearPDEs import LinearPDE, TransportPDE
3  from esys.escript.pdetools import Projector  from esys.escript.pdetools import Projector
4  from esys import finley  from esys import finley
5  import math  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., location_fixed_phi=Data()):  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           """           """
234           self.domain = phi.getDomain()           self.__domain = phi.getDomain()
235           self.__h = Function(self.domain).getSize()           x=self.__domain.getX()
236           self.__location_fixed_phi=location_fixed_phi           diam=0
237             for i in range(self.__domain.getDim()):
238           self.__relative_smoothing_width = relative_smoothing_width               xi=x[i]
239           self.__pde = LinearPDE(self.domain)                       diam+=(inf(xi)-sup(xi))**2
240           self.__pde.setReducedOrderOn()           self.__diam=sqrt(diam)
241           self.__pde.setValue(D=1.0,q=self.__location_fixed_phi)           self.__h = sup(Function(self.__domain).getSize())
242           # self.__pde.setTolerance(1.e-10)           self.__reinit_each=reinit_each
243           #self.__pde.setSolverMethod(solver=LinearPDE.DIRECT)           self.__reinit_max = reinit_max
244                                 self.__smooth = smooth
245           self.__reinitPde = LinearPDE(self.domain,numEquations=1)           self.__phi = phi
          # self.__reinitPde.setSolverMethod(solver=LinearPDE.LUMPING)        
          self.__reinitPde.setReducedOrderOn()  
          self.__pde.setValue(q=self.__location_fixed_phi)  
          self.__reinitPde.setSymmetryOn()  
          self.__reinitPde.setTolerance(1.e-8)  
          # self.__reinitPde.setSolverMethod(preconditioner=LinearPDE.ILU0)  
           
          self.__reinitialization_each=reinitialization_each  
          self.__reinitialization_steps_max = reinitialization_steps_max  
          print self.__reinitialization_steps_max  
246           self.__update_count=0           self.__update_count=0
247             self.velocity = None
248    
249             #==================================================
250           self.__phi = phi           self.__FC=True
251             if self.__FC:
252                self.__fc=TransportPDE(self.__domain,num_equations=1,theta=0.5)
253                # self.__fc.setReducedOrderOn()
254                self.__fc.setValue(M=Scalar(1.,Function(self.__domain)))
255                self.__fc.setInitialSolution(phi+self.__diam)
256             else:
257                self.__fcpde = LinearPDE(self.__domain,numEquations=1, numSolutions=1)
258                self.__fcpde.setSolverMethod(solver=LinearPDE.LUMPING)      
259                self.__fcpde.setReducedOrderOn()
260                self.__fcpde.setSymmetryOn()
261                self.__fcpde.setValue(D=1.)
262    
263             #=======================================
264             self.__reinitFC=False
265             if self.__reinitFC:
266                self.__reinitfc=TransportPDE(self.__domain,num_equations=1,theta=1.0)
267                self.__reinitfc.setValue(M=Scalar(1.,Function(self.__domain)))
268                self.__reinitfc.setTolerance(1.e-5)
269             else:
270                self.__reinitPde = LinearPDE(self.__domain,numEquations=1, numSolutions=1)
271                # self.__reinitPde.setSolverMethod(solver=LinearPDE.LUMPING)      
272                self.__reinitPde.setReducedOrderOn()
273                self.__reinitPde.setSymmetryOn()
274                self.__reinitPde.setValue(D=1.)
275                self.__reinitPde.setTolerance(1.e-2)
276                # self.__reinitPde.setSolverMethod(preconditioner=LinearPDE.ILU0)
277             #=======================================
278           self.__updateInterface()           self.__updateInterface()
279           self.velocity = None           print "phi range:",inf(phi), sup(phi)
280    
281         def setFixedRegion(self,mask,contour=0):
282             q=whereNonPositive(abs(self.__phi-contour)-self.__h)*mask
283             # self.__fc.setValue(q=q)
284             self.__reinitPde.setValue(q=q)
285    
286    
      def getDomain(self):  
          return self.domain  
       
       
287       def __updateInterface(self):       def __updateInterface(self):
288           self.__smoothed_char=self.__makeInterface(self.__relative_smoothing_width)           self.__smoothed_char=self.__makeInterface(self.__phi)
289    
290    
      def __makeInterface(self,smoothing_width=1.):  
          # s=smoothing_width*self.__h*length(grad(self.__phi,self.__h.getFunctionSpace()))  
          s=smoothing_width*self.__h  
          phi_on_h=interpolate(self.__phi,self.__h.getFunctionSpace())          
          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  
291    
      def getCharacteristicFunction(self):  
          return self.__smoothed_char  
292    
293       def update(self,dt):       def update(self,dt):
294           """           """
# Line 67  class LevelSet(object): Line 296  class LevelSet(object):
296    
297           @param dt: time step forward           @param dt: time step forward
298           """           """
299           self.__advectPhi(dt)           if self.__FC:
300           self.__updateInterface()              self.__phi=self.__fc.solve(dt, verbose=False)-self.__diam
301           if self.__update_count%self.__reinitialization_each == 0:           else:
302              self.__reinitialise()              self.__fcpde.setValue(Y = self.__phi-(dt/2.)*inner(self.velocity,grad(self.__phi)))
303              self.__updateInterface()              phi_half = self.__fcpde.getSolution()
304                self.__fcpde.setValue(Y = self.__phi-dt*inner(self.velocity,grad(phi_half)))
305                self.__phi= self.__fcpde.getSolution()
306           self.__update_count += 1           self.__update_count += 1
307             if self.__update_count%self.__reinit_each == 0:
308                self.__phi=self.__reinitialise(self.__phi)
309                if self.__FC:
310                    self.__fc.setInitialSolution(self.__phi+self.__diam)
311             self.__updateInterface()
312    
313         def setTolerance(self,tolerance=1e-3):
314             if self.__FC:
315                self.__fc.setTolerance(tolerance)
316    
317         def __reinitialise(self,phi):
318             print "reintialization started:"
319             s=self.__makeInterface(phi,1.)
320             print "phi range:",inf(phi), sup(phi)
321             g=grad(phi)
322             w=s*g/(length(g)+EPSILON)
323             #=======================================================
324             # # positive part:
325             # phi_p=wherePositive(phi)*phi
326             # self.__reinitfc.setInitialSolution(phi_p)
327             # self.__reinitfc.setValue(C=-wherePositive(s)*w,Y=wherePositive(s)*s,q=whereNonPositive(phi))
328             # dtau=self.__h
329             # print "step size: dt (pos)= ",dtau
330             # print "phi_p range:",inf(phi_p), sup(phi_p)
331             # iter=0
332             # while (iter<=self.__reinit_max):
333             # phi_p=self.__reinitfc.solve(dtau)
334             # print "phi_p range:",inf(phi_p), sup(phi_p)
335             # iter+=1
336             # # negative part:
337             # phi_n=-whereNegative(phi)*phi
338             # self.__reinitfc.setInitialSolution(phi_n)
339             # self.__reinitfc.setValue(C=-whereNegative(s)*w,Y=-whereNegative(s)*s,q=whereNonNegative(phi))
340             # # dtau=self.__reinitfc.getSafeTimeStepSize()
341             # dtau=self.__h
342             # print "step size: dt (neg)= ",dtau
343             # print "phi_n range:",inf(phi_n), sup(phi_n)
344             # iter=0
345             # while (iter<=self.__reinit_max):
346             # phi_n=self.__reinitfc.solve(dtau)
347             # print "phi_n range:",inf(phi_n), sup(phi_n)
348             # iter+=1
349             # phi=phi_p-phi_n
350             # print "phi range:",inf(phi), sup(phi)
351             # print "reintialization completed."
352             # return phi
353    
354             #=======================================================
355             if self.__reinitFC:
356                 self.__reinitfc.setValue(C=-w,Y=s)
357                 self.__reinitfc.setInitialSolution(phi+self.__diam)
358                 dtau=self.__reinitfc.getSafeTimeStepSize()
359                 # dtau = 0.5*inf(Function(self.__domain).getSize())
360             else:
361                 dtau = 0.5*inf(Function(self.__domain).getSize())
362             print "step size: dt = ",dtau,inf(abs(phi.interpolate(Function(self.__domain)))/abs(s-inner(w,grad(phi))))
363             iter =0
364             # self.__reinitPde.setValue(q=whereNegative(abs(phi)-2*self.__h), r=phi)
365             # self.__reinitPde.setValue(r=phi)
366             while (iter<=self.__reinit_max):
367                     phi_old=phi
368                     if self.__reinitFC:
369                       phi = self.__reinitfc.solve(dtau)-self.__diam
370                     else:
371                       self.__reinitPde.setValue(Y = phi+(dtau/2.)*(s-inner(w,grad(phi))))
372                       phi_half = self.__reinitPde.getSolution()
373                       self.__reinitPde.setValue(Y = phi+dtau*(s-inner(w,grad(phi_half))))
374                       # g=grad(phi)
375                       # S=inner(w,grad(phi))
376                       # self.__reinitPde.setValue(Y = phi+dtau*(s-S),X=dtau**2/2*w*(s-S))
377                       phi = self.__reinitPde.getSolution()
378                     change = Lsup(phi-phi_old)/self.__diam
379                     print "phi range:",inf(phi), sup(phi)
380                     print "iteration :", iter, " change:", change
381                     iter +=1
382             print "reintialization completed."
383             return phi
384    
385         def createParameter(self,value_negative=-1.,value_positive=1):
386             out = (value_negative+value_positive)/2. + self.__smoothed_char * ((value_positive-value_negative)/2.)
387             return out
388    
389    
390         #================ things from here onwards are not used nor tested: ==========================================
391        
392    
393         def getCharacteristicFunction(self):
394             return self.__smoothed_char
395    
396    
      def __advectPhi(self,dt):  
         """  
         advects the level set function  
         """  
         fdf=0.01  
         grad_phi = grad(self.__phi)    
         len_grad_phi=length(grad_phi)  
         f=whereNonNegative(len_grad_phi-fdf)  
         s=(1.-f)*len_grad_phi/fdf+f  
           
         # self.__pde.setValue(Y=self.__phi-dt/2.*s*inner(self.velocity,grad_phi),r=self.__phi)    
         # phi_half = self.__pde.getSolution()  
         # self.__pde.setValue(Y=self.__phi-dt*s*inner(self.velocity,grad(phi_half)),r=self.__phi)    
         # self.__phi = self.__pde.getSolution()  
   
         self.__pde.setValue(Y=s*inner(self.velocity,grad_phi))  
         phi_half = self.__phi-dt/2*self.__pde.getSolution()  
         self.__pde.setValue(Y=s*inner(self.velocity,grad(phi_half)))  
         self.__phi =self.__phi-dt*self.__pde.getSolution()  
         print("level set advection done")  
397    
      def getTimeStepSize(self,velocity,safety_factor=0.6):  
          """  
          returns a new dt for a given velocity using the courant coundition  
          """  
          self.velocity=velocity  
          return safety_factor*inf(self.__h/length(interpolate(velocity,self.__h.getFunctionSpace())))  
398    
399       def RK2(self,L):       def RK2(self,L):
400             k0=L(phi)             k0=L(phi)
# Line 121  class LevelSet(object): Line 415  class LevelSet(object):
415             phi=phi+dt*k1             phi=phi+dt*k1
416                        
417                        
418       def __reinitialise(self):       def __reinitialise_old(self):
419          #=============================================          #=============================================
420          f=0.1          f=0.1
421          f2=0.3          f2=0.3
# Line 151  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)
# Line 163  class LevelSet(object): Line 457  class LevelSet(object):
457            # 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))
458            # phi_half=self.__reinitPde.getSolution()            # phi_half=self.__reinitPde.getSolution()
459            # 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))))
460            # self.__reinitPde.setValue(Y =self.__phi, Y_reduced=dtau*s*(1.-inner(w,grad(phi_half,ReducedFunction(self.domain)))))            # self.__reinitPde.setValue(Y =self.__phi, Y_reduced=dtau*s*(1.-inner(w,grad(phi_half,ReducedFunction(self.__domain)))))
461            # self.__reinitPde.setValue(Y =self.__phi+dtau*(s-inner(w,grad_phi)),X=dtau/2*h*(s-inner(w,grad_phi))*w)            # self.__reinitPde.setValue(Y =self.__phi+dtau*(s-inner(w,grad_phi)),X=dtau/2*h*(s-inner(w,grad_phi))*w)
462            # self.__reinitPde.setValue(Y =self.__phi+dtau*s*(1.-len_grad_phi),X=-dtau/2*h*s**2*grad_phi)            # self.__reinitPde.setValue(Y =self.__phi+dtau*s*(1.-len_grad_phi),X=-dtau/2*h*s**2*grad_phi)
463            self.__reinitPde.setValue(Y=self.__phi+dtau*s*(1.-len_grad_phi),X=f*dtau/2*h*abs(s)*(1.-len_grad_phi)*grad_phi/len_grad_phi)            self.__reinitPde.setValue(Y=self.__phi+dtau*s*(1.-len_grad_phi),X=f*dtau/2*h*abs(s)*(1.-len_grad_phi)*grad_phi/len_grad_phi)
464    
465            # self.__reinitPde.setValue(Y=self.__phi+dtau*s*(1.-len_grad_phi),X=f*dtau/2*h*abs(s)*(1.-len_grad_phi)*grad_phi/len_grad_phi)            # self.__reinitPde.setValue(Y=self.__phi+dtau*s*(1.-len_grad_phi),X=f*dtau/2*h*abs(s)*(1.-len_grad_phi)*grad_phi/len_grad_phi)
466            self.__phi, previous = self.__reinitPde.getSolution(), self.__phi            self.__phi, previous = self.__reinitPde.getSolution(), self.__phi
467            self.__updateInterface()            # self.__updateInterface()
468    
469            vol,vol_old=self.getVolumeOfNegativeDomain(),vol            vol,vol_old=self.getVolumeOfNegativeDomain(),vol
470            diff=(vol-vol_old)            diff=(vol-vol_old)
# Line 224  class LevelSet(object): Line 518  class LevelSet(object):
518           """           """
519           return integrate((1.-self.__makeInterface(1.))/2.)           return integrate((1.-self.__makeInterface(1.))/2.)
520    
      def getLevelSetFunction(self):  
          """  
          returns the level set function  
          """  
          return self.__phi  
521                                
522       def getBoundingBoxOfNegativeDomain(self):       def getBoundingBoxOfNegativeDomain(self):
523           """           """
# Line 238  class LevelSet(object): Line 527  class LevelSet(object):
527           mask_phi1=wherePositive(interpolate(self.__phi,fs))           mask_phi1=wherePositive(interpolate(self.__phi,fs))
528           mask_phi2=wherePositive(self.__phi)           mask_phi2=wherePositive(self.__phi)
529           x1=fs.getX()           x1=fs.getX()
530           x2=self.domain.getX()           x2=self.__domain.getX()
531           out=[]           out=[]
532           for i in range(fs.getDim()):           for i in range(fs.getDim()):
533               x1_i=x1[i]               x1_i=x1[i]
# Line 249  class LevelSet(object): Line 538  class LevelSet(object):
538               out.append((min(inf(x1_i+offset1),inf(x2_i+offset2)),max(sup(x1_i-offset1),sup(x2_i-offset2))))               out.append((min(inf(x1_i+offset1),inf(x2_i+offset2)),max(sup(x1_i-offset1),sup(x2_i-offset2))))
539           return tuple(out)           return tuple(out)
540    
      def createParameter(self,value_negative=-1.,value_positive=1):  
          out = (value_negative+value_positive)/2. + self.__smoothed_char * ((value_positive-value_negative)/2.)  
          return out  

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

  ViewVC Help
Powered by ViewVC 1.1.26