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

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

Parent Directory Parent Directory | Revision Log Revision Log


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

  ViewVC Help
Powered by ViewVC 1.1.26