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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1811 - (hide annotations)
Thu Sep 25 23:11:13 2008 UTC (11 years ago) by ksteube
File MIME type: text/x-python
File size: 20945 byte(s)
Copyright updated in all files

1 ksteube 1811
2     ########################################################
3     #
4     # Copyright (c) 2003-2008 by University of Queensland
5     # Earth Systems Science Computational Center (ESSCC)
6     # http://www.uq.edu.au/esscc
7     #
8     # Primary Business: Queensland, Australia
9     # Licensed under the Open Software License version 3.0
10     # http://www.opensource.org/licenses/osl-3.0.php
11     #
12     ########################################################
13    
14     __copyright__="""Copyright (c) 2003-2008 by University of Queensland
15     Earth Systems Science Computational Center (ESSCC)
16     http://www.uq.edu.au/esscc
17     Primary Business: Queensland, Australia"""
18     __license__="""Licensed under the Open Software License version 3.0
19     http://www.opensource.org/licenses/osl-3.0.php"""
20     __url__="http://www.uq.edu.au/esscc/escript-finley"
21    
22 gross 1473 from esys.escript import *
23 gross 1659 from esys.escript.linearPDEs import LinearPDE, TransportPDE
24 gross 1473 from esys.escript.pdetools import Projector
25     from esys import finley
26     import math
27    
28    
29    
30 gross 1788 from esys.escript import *
31     from esys.finley import finley
32     from esys.escript.linearPDEs import LinearPDE
33     from esys.escript.pdetools import Projector
34     import sys
35    
36    
37     class LevelSet:
38     """
39     The level set method tracking an interface defined by the zero contour of the level set function phi.
40     it is assumed the phi(x)<0 defines the volume of interest.
41    
42     """
43     def __init__(self,domain,phi,reinit_max=10,reinit_each=2,smooth=2.):
44     """
45     set up the level set method
46    
47     @param domain: the domain where the level set is used
48     @param phi: the initial level set function
49     @param reinit_max: maximum number of reinitalization steps
50     @param reinit_each: phi is reinitialized every reinit_each step
51     @param smooth: smoothing width
52     """
53     self.__domain = domain
54     self.__phi = phi
55     self.__reinit_max = reinit_max
56     self.__reinit_each = reinit_each
57     self.__PDE = LinearPDE(domain)
58     self.__PDE.setReducedOrderOn()
59     self.__PDE.setValue(D=1.0)
60     self.__PDE.setSolverMethod(solver=LinearPDE.PCG)
61     self.__reinitPDE = LinearPDE(domain, numEquations=1)
62     self.__reinitPDE.setReducedOrderOn()
63     self.__reinitPDE.setValue(D=1.0)
64     self.__reinitPDE.setSolverMethod(solver=LinearPDE.LUMPING)
65     self.__h = inf(domain.getSize())
66     self.__smooth = smooth
67     self.__n_step=0
68    
69     def __advect(self, velocity, dt):
70     """
71     advects the level set function in the presense of a velocity field.
72    
73     This implementation uses the 2-step Taylor-Galerkin method
74     @param velocity: velocity field
75     @param dt: dime increment
76     @return: the advected level set function
77     """
78     Y = self.__phi-dt/2.*inner(velocity,grad(self.__phi))
79     self.__PDE.setValue(Y=Y)
80     phi_half = self.__PDE.getSolution()
81     Y = self.__phi-dt*inner(velocity,grad(phi_half))
82     self.__PDE.setValue(Y=Y)
83     phi = self.__PDE.getSolution()
84     print "LevelSet: Advection step done"
85     return phi
86    
87     def __reinitialise(self):
88     """
89     reinializes the level set
90    
91     It solves the
92    
93     @return: reinitalized level set
94     """
95     phi=self.__phi
96     s = sign(phi.interpolate(Function(self.__domain)))
97     g=grad(phi)
98     w = s*g/(length(g)+1e-10)
99     dtau = 0.5*self.__h
100     iter =0
101     mask = whereNegative(abs(phi)-1.*self.__h)
102     self.__reinitPDE.setValue(q=mask, r=phi)
103     g=grad(phi)
104     while (iter<=self.__reinit_max):
105     Y = phi+(dtau/2.)*(s-inner(w,g))
106     self.__reinitPDE.setValue(Y = Y)
107     phi_half = self.__reinitPDE.getSolution()
108     Y = phi+dtau*(s-inner(w,grad(phi_half)))
109     self.__reinitPDE.setValue(Y = Y)
110     phi = self.__reinitPDE.getSolution()
111     g=grad(phi)
112     error = Lsup(length(g)*whereNegative(abs(phi.interpolate(g.getFunctionSpace()))-3.0*self.__h))
113     print "LevelSet:reinitialization: iteration :", iter, " error:", error
114     iter +=1
115     return phi
116    
117     def getTimeStepSize(self,velocity):
118     """
119     returns a new dt for a given velocity using the courant coundition
120    
121     @param velocity: velocity field
122     """
123     self.__velocity=velocity
124     dt=0.5*self.__h/sup(length(velocity))
125     return dt
126    
127     def update(self,dt):
128     """
129     sets a new velocity and updates the level set fuction
130    
131     @param dt: time step forward
132     """
133     self.__phi=self.__advect(self.__velocity, dt)
134     self.__n_step+=1
135     if self.__n_step%self.__reinit_each ==0: self.__phi = self.__reinitialise()
136     return self.__phi
137    
138     def update_phi(self, velocity, dt):
139     """
140     updates phi under the presense of a velocity field
141    
142     If dt is small this call is equivalent to call
143    
144     dt=LevelSet.getTimeStepSize(velocity)
145     phi=LevelSet.update(dt)
146    
147     otherwise substepping is used.
148     @param velocity: velocity field
149     @param dt: time step forward
150     """
151     dt2=self.getTimeStepSize(velocity)
152     n=math.ceil(dt/dt2)
153     dt_new=dt/n
154     for i in range(n):
155     phi=self.update(dt_new)
156     t+=dt_new
157     return phi
158    
159    
160     def getVolume(self):
161     """
162     return the volume of the phi(x)<0 region
163     """
164     return integrate(whereNegative(self.__phi.interpolate(Function(self.__domain))))
165    
166     def getSurface(self,rel_width_factor=0.5):
167     """
168     return a mask for phi(x)=1 region
169    
170     @param rel_width_factor: relative wideth of region around zero contour.
171     """
172     return whereNegative(abs(self.__phi)-rel_width_factor*self.__h)
173    
174     def getH(self):
175     """
176     returns mesh size
177     """
178     return self.__h
179    
180     def getDomain(self):
181     """
182     returns domain
183     """
184     return self.__domain
185    
186     def getLevelSetFunction(self):
187     """
188     returns the level set function
189     """
190     return self.__phi
191    
192     def update_parameter_sharp(self, param_neg=-1, param_pos=1, phi=None):
193     """
194     creates a function whith param_neg where phi<0 and param_pos where phi>0 (no smoothing)
195    
196     @param param_neg: value of parameter on the negative side (phi<0)
197     @param param_pos: value of parameter on the positve side (phi>0)
198     @param phi: level set funtion to be used. if not present the current level set is used.
199     """
200     mask_neg = whereNegative(self.__phi)
201     mask_pos = whereNonNegative(self.__phi)
202     param = param_pos*mask_pos + param_neg*mask_neg
203     return param
204    
205     def update_parameter(self, param_neg=-1, param_pos=1, phi=None, smoothing_width=None):
206     """
207     creates a smoothed function whith param_neg where phi<0 and param_pos where phi>0 which is smoothed over a length
208     smoothing_width accross the interface
209    
210     @param smoothing_width: width of the smoothing zone relative to mesh size. If not present the initial value of C{smooth} is used.
211     """
212     if smoothing_width==None: smoothing_width = self.__smooth
213     if phi==None: phi=self.__phi
214     s=self.__makeInterface(phi,smoothing_width)
215     return ((param_pos-param_neg)*s+param_pos+param_neg)/2
216    
217     def __makeInterface(self,phi,smoothing_width):
218     """
219     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
220     and 1 where the level set is 1
221     """
222     s=smoothing_width*self.__h
223     phi_on_h=interpolate(phi,Function(self.__domain))
224     mask_neg = whereNonNegative(-s-phi_on_h)
225     mask_pos = whereNonNegative(phi_on_h-s)
226     mask_interface = 1.-mask_neg-mask_pos
227     interface=phi_on_h/s
228     return - mask_neg + mask_pos + mask_interface * interface
229    
230     def makeCharacteristicFunction(self, contour=0, phi=None, positiveSide=True, smoothing_width=None):
231     """
232     makes a smooth charateristic function of the region phi(x)>contour if positiveSide and phi(x)<contour otherwise.
233    
234     @param phi: level set funtion to be used. if not present the current level set is used.
235     @param smoothing_width: width of the smoothing zone relative to mesh size. If not present the initial value of C{smooth} is used.
236     """
237     if phi==None: phi=self.__phi
238     if smoothing_width == None: smoothing_width=self.__smooth
239     s=self.__makeInterface(phi=phi-contour,smoothing_width=smoothing_width)
240     if positiveSide:
241     return (1+s)/2
242     else:
243     return (1-s)/2
244    
245     def setTolerance(self,tolerance=1e-3):
246     self.__PDE.setTolerance(tolerance)
247     self.__reinitPDE.setTolerance(tolerance)
248    
249    
250     class LevelSet2(object):
251     def __init__(self,phi,reinit_max=10,reinit_each=2,smooth=2.):
252 gross 1473 """
253     initialize model
254     """
255 gross 1659 self.__domain = phi.getDomain()
256     x=self.__domain.getX()
257     diam=0
258     for i in range(self.__domain.getDim()):
259     xi=x[i]
260     diam+=(inf(xi)-sup(xi))**2
261     self.__diam=sqrt(diam)
262 gross 1661 self.__h = sup(Function(self.__domain).getSize())
263 gross 1788 self.__reinit_each=reinit_each
264     self.__reinit_max = reinit_max
265     self.__smooth = smooth
266 gross 1659 self.__phi = phi
267     self.__update_count=0
268     self.velocity = None
269 gross 1473
270 gross 1659 #==================================================
271 gross 1673 self.__FC=True
272     if self.__FC:
273     self.__fc=TransportPDE(self.__domain,num_equations=1,theta=0.5)
274     # self.__fc.setReducedOrderOn()
275     self.__fc.setValue(M=Scalar(1.,Function(self.__domain)))
276     self.__fc.setInitialSolution(phi+self.__diam)
277     else:
278     self.__fcpde = LinearPDE(self.__domain,numEquations=1, numSolutions=1)
279     self.__fcpde.setSolverMethod(solver=LinearPDE.LUMPING)
280     self.__fcpde.setReducedOrderOn()
281     self.__fcpde.setSymmetryOn()
282     self.__fcpde.setValue(D=1.)
283 gross 1659
284     #=======================================
285 gross 1673 self.__reinitFC=False
286     if self.__reinitFC:
287 gross 1716 self.__reinitfc=TransportPDE(self.__domain,num_equations=1,theta=1.0)
288 gross 1673 self.__reinitfc.setValue(M=Scalar(1.,Function(self.__domain)))
289     self.__reinitfc.setTolerance(1.e-5)
290     else:
291     self.__reinitPde = LinearPDE(self.__domain,numEquations=1, numSolutions=1)
292     # self.__reinitPde.setSolverMethod(solver=LinearPDE.LUMPING)
293     self.__reinitPde.setReducedOrderOn()
294     self.__reinitPde.setSymmetryOn()
295     self.__reinitPde.setValue(D=1.)
296 gross 1716 self.__reinitPde.setTolerance(1.e-2)
297 gross 1673 # self.__reinitPde.setSolverMethod(preconditioner=LinearPDE.ILU0)
298 gross 1659 #=======================================
299 gross 1661 self.__updateInterface()
300     print "phi range:",inf(phi), sup(phi)
301 gross 1473
302 gross 1661 def setFixedRegion(self,mask,contour=0):
303     q=whereNonPositive(abs(self.__phi-contour)-self.__h)*mask
304     # self.__fc.setValue(q=q)
305     self.__reinitPde.setValue(q=q)
306 gross 1473
307    
308 gross 1661 def __updateInterface(self):
309     self.__smoothed_char=self.__makeInterface(self.__phi)
310    
311 gross 1473
312 gross 1661
313    
314 gross 1473 def update(self,dt):
315     """
316     sets a new velocity and updates the level set fuction
317    
318     @param dt: time step forward
319     """
320 gross 1673 if self.__FC:
321     self.__phi=self.__fc.solve(dt, verbose=False)-self.__diam
322     else:
323     self.__fcpde.setValue(Y = self.__phi-(dt/2.)*inner(self.velocity,grad(self.__phi)))
324     phi_half = self.__fcpde.getSolution()
325     self.__fcpde.setValue(Y = self.__phi-dt*inner(self.velocity,grad(phi_half)))
326     self.__phi= self.__fcpde.getSolution()
327 gross 1659 self.__update_count += 1
328 gross 1788 if self.__update_count%self.__reinit_each == 0:
329 gross 1661 self.__phi=self.__reinitialise(self.__phi)
330 gross 1673 if self.__FC:
331     self.__fc.setInitialSolution(self.__phi+self.__diam)
332 gross 1661 self.__updateInterface()
333 gross 1473
334 gross 1659 def setTolerance(self,tolerance=1e-3):
335 gross 1673 if self.__FC:
336     self.__fc.setTolerance(tolerance)
337 gross 1473
338 gross 1659 def __reinitialise(self,phi):
339     print "reintialization started:"
340 gross 1716 s=self.__makeInterface(phi,1.)
341     print "phi range:",inf(phi), sup(phi)
342 gross 1659 g=grad(phi)
343 gross 1716 w=s*g/(length(g)+EPSILON)
344     #=======================================================
345     # # positive part:
346     # phi_p=wherePositive(phi)*phi
347     # self.__reinitfc.setInitialSolution(phi_p)
348     # self.__reinitfc.setValue(C=-wherePositive(s)*w,Y=wherePositive(s)*s,q=whereNonPositive(phi))
349     # dtau=self.__h
350     # print "step size: dt (pos)= ",dtau
351     # print "phi_p range:",inf(phi_p), sup(phi_p)
352     # iter=0
353 gross 1788 # while (iter<=self.__reinit_max):
354 gross 1716 # phi_p=self.__reinitfc.solve(dtau)
355     # print "phi_p range:",inf(phi_p), sup(phi_p)
356     # iter+=1
357     # # negative part:
358     # phi_n=-whereNegative(phi)*phi
359     # self.__reinitfc.setInitialSolution(phi_n)
360     # self.__reinitfc.setValue(C=-whereNegative(s)*w,Y=-whereNegative(s)*s,q=whereNonNegative(phi))
361     # # dtau=self.__reinitfc.getSafeTimeStepSize()
362     # dtau=self.__h
363     # print "step size: dt (neg)= ",dtau
364     # print "phi_n range:",inf(phi_n), sup(phi_n)
365     # iter=0
366 gross 1788 # while (iter<=self.__reinit_max):
367 gross 1716 # phi_n=self.__reinitfc.solve(dtau)
368     # print "phi_n range:",inf(phi_n), sup(phi_n)
369     # iter+=1
370     # phi=phi_p-phi_n
371     # print "phi range:",inf(phi), sup(phi)
372     # print "reintialization completed."
373     # return phi
374    
375     #=======================================================
376 gross 1673 if self.__reinitFC:
377     self.__reinitfc.setValue(C=-w,Y=s)
378     self.__reinitfc.setInitialSolution(phi+self.__diam)
379 gross 1716 dtau=self.__reinitfc.getSafeTimeStepSize()
380     # dtau = 0.5*inf(Function(self.__domain).getSize())
381 gross 1673 else:
382     dtau = 0.5*inf(Function(self.__domain).getSize())
383     print "step size: dt = ",dtau,inf(abs(phi.interpolate(Function(self.__domain)))/abs(s-inner(w,grad(phi))))
384 gross 1659 iter =0
385 gross 1673 # self.__reinitPde.setValue(q=whereNegative(abs(phi)-2*self.__h), r=phi)
386 gross 1661 # self.__reinitPde.setValue(r=phi)
387 gross 1788 while (iter<=self.__reinit_max):
388 gross 1659 phi_old=phi
389 gross 1673 if self.__reinitFC:
390     phi = self.__reinitfc.solve(dtau)-self.__diam
391     else:
392     self.__reinitPde.setValue(Y = phi+(dtau/2.)*(s-inner(w,grad(phi))))
393     phi_half = self.__reinitPde.getSolution()
394     self.__reinitPde.setValue(Y = phi+dtau*(s-inner(w,grad(phi_half))))
395     # g=grad(phi)
396     # S=inner(w,grad(phi))
397     # self.__reinitPde.setValue(Y = phi+dtau*(s-S),X=dtau**2/2*w*(s-S))
398     phi = self.__reinitPde.getSolution()
399 gross 1659 change = Lsup(phi-phi_old)/self.__diam
400     print "phi range:",inf(phi), sup(phi)
401 gross 1716 print "iteration :", iter, " change:", change
402 gross 1659 iter +=1
403     print "reintialization completed."
404     return phi
405 gross 1473
406 gross 1661 def createParameter(self,value_negative=-1.,value_positive=1):
407     out = (value_negative+value_positive)/2. + self.__smoothed_char * ((value_positive-value_negative)/2.)
408     return out
409 gross 1473
410 gross 1661
411 gross 1659 #================ things from here onwards are not used nor tested: ==========================================
412    
413    
414     def getCharacteristicFunction(self):
415     return self.__smoothed_char
416    
417    
418    
419    
420 gross 1473 def RK2(self,L):
421     k0=L(phi)
422     phi_1=phi+dt*k0
423     k1=L(phi_1)
424     phi=phi+dt/2*(k0+k1)
425     def RK3(self,L):
426     k0=L(phi)
427     phi_1=phi+dt*L(phi)
428     k1=L(phi_1)
429     phi_2=phi+dt/4*(k0+k1)
430     k2=L(phi_2)
431     phi=phi+dt/6*(k0+4*k2+K1)
432     def TG(self,L):
433     k0=L(phi)
434     phi_1=phi+dt/2*k0
435     k1=L(phi_1)
436     phi=phi+dt*k1
437    
438    
439 gross 1659 def __reinitialise_old(self):
440 gross 1473 #=============================================
441     f=0.1
442     f2=0.3
443     h=self.__h
444     fs=h.getFunctionSpace()
445     vol=self.getVolumeOfNegativeDomain()
446     self.__reinitPde.setValue(D=1.0)
447     #============
448     grad_phi=grad(self.__phi,fs)
449     len_grad_phi=length(grad_phi)
450     w=grad_phi/len_grad_phi
451    
452     # s=wherePositive(self. __makeInterface(2.))-whereNegative(self. __makeInterface(2.))
453     s=self. __makeInterface(2.)
454     # s2=whereNegative(len_grad_phi-f2)*(1-(1-len_grad_phi/f2)**2)+(1-whereNegative(len_grad_phi-f2))
455     s2=1.
456     s=s*s2
457     dtau=f*inf(h/abs(s))
458    
459     #========================
460     diff=1.
461     d=1.
462     c =0
463     #==============================================
464     TVD=integrate(length(grad_phi))
465     print "initial range ",inf(self.__phi),sup(self.__phi),"error:",Lsup(1.-len_grad_phi),"volume =",vol,TVD
466     # saveVTK("test.%s.xml"%c,l=length(grad(self.__phi,fs))-1,s=s,phi=self.__phi)
467    
468     dtau=f*inf(h/abs(s))
469 gross 1788 while c < self.__reinit_max: # and abs(diff) >= 0.01:
470 gross 1473 #
471     grad_phi=grad(self.__phi,fs)
472     len_grad_phi=length(grad_phi)
473     #
474     # s=self.__makeInterface(1.)
475     # s2=whereNegative(len_grad_phi-f2)*(1-(1-len_grad_phi/f2)**2)+(1-whereNegative(len_grad_phi-f2))
476     # s=s*s2
477     # phi_on_h=interpolate(self.__phi,fs)
478     # self.__reinitPde.setValue(Y =self.__phi+dtau/2*s*(1.-len_grad_phi))
479     # phi_half=self.__reinitPde.getSolution()
480     # self.__reinitPde.setValue(Y =self.__phi+dtau*s*(1.-length(grad(phi_half,fs))))
481 gross 1659 # self.__reinitPde.setValue(Y =self.__phi, Y_reduced=dtau*s*(1.-inner(w,grad(phi_half,ReducedFunction(self.__domain)))))
482 gross 1473 # self.__reinitPde.setValue(Y =self.__phi+dtau*(s-inner(w,grad_phi)),X=dtau/2*h*(s-inner(w,grad_phi))*w)
483     # self.__reinitPde.setValue(Y =self.__phi+dtau*s*(1.-len_grad_phi),X=-dtau/2*h*s**2*grad_phi)
484     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)
485    
486     # 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)
487     self.__phi, previous = self.__reinitPde.getSolution(), self.__phi
488 gross 1661 # self.__updateInterface()
489 gross 1473
490     vol,vol_old=self.getVolumeOfNegativeDomain(),vol
491     diff=(vol-vol_old)
492     r=Lsup(length(grad(self.__phi))-1.)
493     TVD=integrate(length(grad(self.__phi,fs)))
494     print "iteration :", c, "range ",inf(self.__phi),sup(self.__phi),"error :",r,"volume change:",diff,TVD
495     # saveVTK("test.%s.xml"%(c+1),l=length(grad(self.__phi,fs)),s=s,phi=self.__phi,v=grad(self.__phi,fs))
496     c += 1
497     return
498     #==============================================
499     f2=0.7
500     h=self.__h
501     fs=h.getFunctionSpace()
502     vol=self.getVolumeOfNegativeDomain()
503     s=abs(self. __makeInterface(2.))
504     grad_phi=grad(self.__phi,fs)
505     len_grad_phi=length(grad_phi)
506    
507     #----
508     # aphi_on_h=abs(interpolate(self.__phi,self.__h.getFunctionSpace()))
509     # q=2*self.__h
510     # mask_neg = whereNegative(aphi_on_h-q)
511     # mask_pos = wherePositive(aphi_on_h-q-self.__h)
512     # mask_interface = 1.-mask_neg-mask_pos
513     # s2=mask_neg + mask_interface * (q+self.__h-aphi_on_h)/self.__h
514     #----
515     m=whereNonPositive(len_grad_phi-f2)
516     s2=m*len_grad_phi/f2+(1.-m)
517     #----
518     # s2=1.
519     #----
520     self.__reinitPde.setValue(D=(1-s)/h**2,Y=(1-s)/h**2*self.__phi)
521     self.__reinitPde.setValue(A=(s2*len_grad_phi+(1.-s2))*kronecker(3),X=grad_phi)
522     # self.__reinitPde.setValue(A=(len_grad_phi-1)/len_grad_phi*kronecker(3))
523     # self.__reinitPde.setValue(A=kronecker(3), X=grad_phi/len_grad_phi)
524     self.__phi = self.__reinitPde.getSolution()
525     r=Lsup(length(grad(self.__phi))-1.)
526     vol,vol_old=self.getVolumeOfNegativeDomain(),vol
527     diff=(vol-vol_old)/vol
528     print "iteration :", inf(self.__phi),sup(self.__phi),r,diff
529     # saveVTK("test.%s.xml"%0,l=length(grad(self.__phi,fs)),s=s,phi=self.__phi,v=grad(self.__phi,fs),s2=s2)
530     return
531     #=============================================
532    
533     #============
534    
535    
536     def getVolumeOfNegativeDomain(self):
537     """
538     return the current volume of domain with phi<0.
539     """
540     return integrate((1.-self.__makeInterface(1.))/2.)
541    
542    
543     def getBoundingBoxOfNegativeDomain(self):
544     """
545     get the height of the region with phi<0
546     """
547     fs=self.__h.getFunctionSpace()
548     mask_phi1=wherePositive(interpolate(self.__phi,fs))
549     mask_phi2=wherePositive(self.__phi)
550     x1=fs.getX()
551 gross 1659 x2=self.__domain.getX()
552 gross 1473 out=[]
553     for i in range(fs.getDim()):
554     x1_i=x1[i]
555     x2_i=x2[i]
556     d=2*(sup(x2_i)-inf(x2_i))
557     offset1=d*mask_phi1
558     offset2=d*mask_phi2
559     out.append((min(inf(x1_i+offset1),inf(x2_i+offset2)),max(sup(x1_i-offset1),sup(x2_i-offset2))))
560     return tuple(out)
561    

  ViewVC Help
Powered by ViewVC 1.1.26