/[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 1788 - (show annotations)
Mon Sep 15 02:20:20 2008 UTC (11 years, 1 month ago) by gross
File MIME type: text/x-python
File size: 20180 byte(s)
levelset in the way laurent did it.
1 from esys.escript import *
2 from esys.escript.linearPDEs import LinearPDE, TransportPDE
3 from esys.escript.pdetools import Projector
4 from esys import finley
5 import math
6
7
8
9 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 """
232 initialize model
233 """
234 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 self.__h = sup(Function(self.__domain).getSize())
242 self.__reinit_each=reinit_each
243 self.__reinit_max = reinit_max
244 self.__smooth = smooth
245 self.__phi = phi
246 self.__update_count=0
247 self.velocity = None
248
249 #==================================================
250 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()
279 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
287 def __updateInterface(self):
288 self.__smoothed_char=self.__makeInterface(self.__phi)
289
290
291
292
293 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 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 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
397
398
399 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 def __reinitialise_old(self):
419 #=============================================
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 while c < self.__reinit_max: # and abs(diff) >= 0.01:
449 #
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 # 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)
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 # self.__updateInterface()
468
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 x2=self.__domain.getX()
531 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