/[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 2344 - (show annotations)
Mon Mar 30 02:13:58 2009 UTC (10 years, 5 months ago) by jfenwick
File MIME type: text/x-python
File size: 21146 byte(s)
Change __url__ to launchpad site

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

  ViewVC Help
Powered by ViewVC 1.1.26