/[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 1788 - (hide annotations)
Mon Sep 15 02:20:20 2008 UTC (11 years, 5 months ago) by gross
File MIME type: text/x-python
File size: 20180 byte(s)
levelset in the way laurent did it.
1 gross 1473 from esys.escript import *
2 gross 1659 from esys.escript.linearPDEs import LinearPDE, TransportPDE
3 gross 1473 from esys.escript.pdetools import Projector
4     from esys import finley
5     import math
6    
7    
8    
9 gross 1788 from esys.escript import *
10     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 gross 1473 """
232     initialize model
233     """
234 gross 1659 self.__domain = phi.getDomain()
235     x=self.__domain.getX()
236     diam=0
237     for i in range(self.__domain.getDim()):
238     xi=x[i]
239     diam+=(inf(xi)-sup(xi))**2
240     self.__diam=sqrt(diam)
241 gross 1661 self.__h = sup(Function(self.__domain).getSize())
242 gross 1788 self.__reinit_each=reinit_each
243     self.__reinit_max = reinit_max
244     self.__smooth = smooth
245 gross 1659 self.__phi = phi
246     self.__update_count=0
247     self.velocity = None
248 gross 1473
249 gross 1659 #==================================================
250 gross 1673 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 gross 1659
263     #=======================================
264 gross 1673 self.__reinitFC=False
265     if self.__reinitFC:
266 gross 1716 self.__reinitfc=TransportPDE(self.__domain,num_equations=1,theta=1.0)
267 gross 1673 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 gross 1716 self.__reinitPde.setTolerance(1.e-2)
276 gross 1673 # self.__reinitPde.setSolverMethod(preconditioner=LinearPDE.ILU0)
277 gross 1659 #=======================================
278 gross 1661 self.__updateInterface()
279     print "phi range:",inf(phi), sup(phi)
280 gross 1473
281 gross 1661 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 gross 1473
286    
287 gross 1661 def __updateInterface(self):
288     self.__smoothed_char=self.__makeInterface(self.__phi)
289    
290 gross 1473
291 gross 1661
292    
293 gross 1473 def update(self,dt):
294     """
295     sets a new velocity and updates the level set fuction
296    
297     @param dt: time step forward
298     """
299 gross 1673 if self.__FC:
300     self.__phi=self.__fc.solve(dt, verbose=False)-self.__diam
301     else:
302     self.__fcpde.setValue(Y = self.__phi-(dt/2.)*inner(self.velocity,grad(self.__phi)))
303     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 gross 1659 self.__update_count += 1
307 gross 1788 if self.__update_count%self.__reinit_each == 0:
308 gross 1661 self.__phi=self.__reinitialise(self.__phi)
309 gross 1673 if self.__FC:
310     self.__fc.setInitialSolution(self.__phi+self.__diam)
311 gross 1661 self.__updateInterface()
312 gross 1473
313 gross 1659 def setTolerance(self,tolerance=1e-3):
314 gross 1673 if self.__FC:
315     self.__fc.setTolerance(tolerance)
316 gross 1473
317 gross 1659 def __reinitialise(self,phi):
318     print "reintialization started:"
319 gross 1716 s=self.__makeInterface(phi,1.)
320     print "phi range:",inf(phi), sup(phi)
321 gross 1659 g=grad(phi)
322 gross 1716 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 gross 1788 # while (iter<=self.__reinit_max):
333 gross 1716 # 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 gross 1788 # while (iter<=self.__reinit_max):
346 gross 1716 # 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 gross 1673 if self.__reinitFC:
356     self.__reinitfc.setValue(C=-w,Y=s)
357     self.__reinitfc.setInitialSolution(phi+self.__diam)
358 gross 1716 dtau=self.__reinitfc.getSafeTimeStepSize()
359     # dtau = 0.5*inf(Function(self.__domain).getSize())
360 gross 1673 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 gross 1659 iter =0
364 gross 1673 # self.__reinitPde.setValue(q=whereNegative(abs(phi)-2*self.__h), r=phi)
365 gross 1661 # self.__reinitPde.setValue(r=phi)
366 gross 1788 while (iter<=self.__reinit_max):
367 gross 1659 phi_old=phi
368 gross 1673 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 gross 1659 change = Lsup(phi-phi_old)/self.__diam
379     print "phi range:",inf(phi), sup(phi)
380 gross 1716 print "iteration :", iter, " change:", change
381 gross 1659 iter +=1
382     print "reintialization completed."
383     return phi
384 gross 1473
385 gross 1661 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 gross 1473
389 gross 1661
390 gross 1659 #================ things from here onwards are not used nor tested: ==========================================
391    
392    
393     def getCharacteristicFunction(self):
394     return self.__smoothed_char
395    
396    
397    
398    
399 gross 1473 def RK2(self,L):
400     k0=L(phi)
401     phi_1=phi+dt*k0
402     k1=L(phi_1)
403     phi=phi+dt/2*(k0+k1)
404     def RK3(self,L):
405     k0=L(phi)
406     phi_1=phi+dt*L(phi)
407     k1=L(phi_1)
408     phi_2=phi+dt/4*(k0+k1)
409     k2=L(phi_2)
410     phi=phi+dt/6*(k0+4*k2+K1)
411     def TG(self,L):
412     k0=L(phi)
413     phi_1=phi+dt/2*k0
414     k1=L(phi_1)
415     phi=phi+dt*k1
416    
417    
418 gross 1659 def __reinitialise_old(self):
419 gross 1473 #=============================================
420     f=0.1
421     f2=0.3
422     h=self.__h
423     fs=h.getFunctionSpace()
424     vol=self.getVolumeOfNegativeDomain()
425     self.__reinitPde.setValue(D=1.0)
426     #============
427     grad_phi=grad(self.__phi,fs)
428     len_grad_phi=length(grad_phi)
429     w=grad_phi/len_grad_phi
430    
431     # s=wherePositive(self. __makeInterface(2.))-whereNegative(self. __makeInterface(2.))
432     s=self. __makeInterface(2.)
433     # s2=whereNegative(len_grad_phi-f2)*(1-(1-len_grad_phi/f2)**2)+(1-whereNegative(len_grad_phi-f2))
434     s2=1.
435     s=s*s2
436     dtau=f*inf(h/abs(s))
437    
438     #========================
439     diff=1.
440     d=1.
441     c =0
442     #==============================================
443     TVD=integrate(length(grad_phi))
444     print "initial range ",inf(self.__phi),sup(self.__phi),"error:",Lsup(1.-len_grad_phi),"volume =",vol,TVD
445     # 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))
448 gross 1788 while c < self.__reinit_max: # and abs(diff) >= 0.01:
449 gross 1473 #
450     grad_phi=grad(self.__phi,fs)
451     len_grad_phi=length(grad_phi)
452     #
453     # s=self.__makeInterface(1.)
454     # s2=whereNegative(len_grad_phi-f2)*(1-(1-len_grad_phi/f2)**2)+(1-whereNegative(len_grad_phi-f2))
455     # s=s*s2
456     # phi_on_h=interpolate(self.__phi,fs)
457     # self.__reinitPde.setValue(Y =self.__phi+dtau/2*s*(1.-len_grad_phi))
458     # phi_half=self.__reinitPde.getSolution()
459     # self.__reinitPde.setValue(Y =self.__phi+dtau*s*(1.-length(grad(phi_half,fs))))
460 gross 1659 # self.__reinitPde.setValue(Y =self.__phi, Y_reduced=dtau*s*(1.-inner(w,grad(phi_half,ReducedFunction(self.__domain)))))
461 gross 1473 # 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)
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)
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)
466     self.__phi, previous = self.__reinitPde.getSolution(), self.__phi
467 gross 1661 # self.__updateInterface()
468 gross 1473
469     vol,vol_old=self.getVolumeOfNegativeDomain(),vol
470     diff=(vol-vol_old)
471     r=Lsup(length(grad(self.__phi))-1.)
472     TVD=integrate(length(grad(self.__phi,fs)))
473     print "iteration :", c, "range ",inf(self.__phi),sup(self.__phi),"error :",r,"volume change:",diff,TVD
474     # saveVTK("test.%s.xml"%(c+1),l=length(grad(self.__phi,fs)),s=s,phi=self.__phi,v=grad(self.__phi,fs))
475     c += 1
476     return
477     #==============================================
478     f2=0.7
479     h=self.__h
480     fs=h.getFunctionSpace()
481     vol=self.getVolumeOfNegativeDomain()
482     s=abs(self. __makeInterface(2.))
483     grad_phi=grad(self.__phi,fs)
484     len_grad_phi=length(grad_phi)
485    
486     #----
487     # aphi_on_h=abs(interpolate(self.__phi,self.__h.getFunctionSpace()))
488     # q=2*self.__h
489     # mask_neg = whereNegative(aphi_on_h-q)
490     # mask_pos = wherePositive(aphi_on_h-q-self.__h)
491     # mask_interface = 1.-mask_neg-mask_pos
492     # s2=mask_neg + mask_interface * (q+self.__h-aphi_on_h)/self.__h
493     #----
494     m=whereNonPositive(len_grad_phi-f2)
495     s2=m*len_grad_phi/f2+(1.-m)
496     #----
497     # s2=1.
498     #----
499     self.__reinitPde.setValue(D=(1-s)/h**2,Y=(1-s)/h**2*self.__phi)
500     self.__reinitPde.setValue(A=(s2*len_grad_phi+(1.-s2))*kronecker(3),X=grad_phi)
501     # self.__reinitPde.setValue(A=(len_grad_phi-1)/len_grad_phi*kronecker(3))
502     # self.__reinitPde.setValue(A=kronecker(3), X=grad_phi/len_grad_phi)
503     self.__phi = self.__reinitPde.getSolution()
504     r=Lsup(length(grad(self.__phi))-1.)
505     vol,vol_old=self.getVolumeOfNegativeDomain(),vol
506     diff=(vol-vol_old)/vol
507     print "iteration :", inf(self.__phi),sup(self.__phi),r,diff
508     # saveVTK("test.%s.xml"%0,l=length(grad(self.__phi,fs)),s=s,phi=self.__phi,v=grad(self.__phi,fs),s2=s2)
509     return
510     #=============================================
511    
512     #============
513    
514    
515     def getVolumeOfNegativeDomain(self):
516     """
517     return the current volume of domain with phi<0.
518     """
519     return integrate((1.-self.__makeInterface(1.))/2.)
520    
521    
522     def getBoundingBoxOfNegativeDomain(self):
523     """
524     get the height of the region with phi<0
525     """
526     fs=self.__h.getFunctionSpace()
527     mask_phi1=wherePositive(interpolate(self.__phi,fs))
528     mask_phi2=wherePositive(self.__phi)
529     x1=fs.getX()
530 gross 1659 x2=self.__domain.getX()
531 gross 1473 out=[]
532     for i in range(fs.getDim()):
533     x1_i=x1[i]
534     x2_i=x2[i]
535     d=2*(sup(x2_i)-inf(x2_i))
536     offset1=d*mask_phi1
537     offset2=d*mask_phi2
538     out.append((min(inf(x1_i+offset1),inf(x2_i+offset2)),max(sup(x1_i-offset1),sup(x2_i-offset2))))
539     return tuple(out)
540    

  ViewVC Help
Powered by ViewVC 1.1.26