/[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 1811 - (show 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
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 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 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 """
253 initialize model
254 """
255 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 self.__h = sup(Function(self.__domain).getSize())
263 self.__reinit_each=reinit_each
264 self.__reinit_max = reinit_max
265 self.__smooth = smooth
266 self.__phi = phi
267 self.__update_count=0
268 self.velocity = None
269
270 #==================================================
271 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
284 #=======================================
285 self.__reinitFC=False
286 if self.__reinitFC:
287 self.__reinitfc=TransportPDE(self.__domain,num_equations=1,theta=1.0)
288 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 self.__reinitPde.setTolerance(1.e-2)
297 # self.__reinitPde.setSolverMethod(preconditioner=LinearPDE.ILU0)
298 #=======================================
299 self.__updateInterface()
300 print "phi range:",inf(phi), sup(phi)
301
302 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
307
308 def __updateInterface(self):
309 self.__smoothed_char=self.__makeInterface(self.__phi)
310
311
312
313
314 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 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 self.__update_count += 1
328 if self.__update_count%self.__reinit_each == 0:
329 self.__phi=self.__reinitialise(self.__phi)
330 if self.__FC:
331 self.__fc.setInitialSolution(self.__phi+self.__diam)
332 self.__updateInterface()
333
334 def setTolerance(self,tolerance=1e-3):
335 if self.__FC:
336 self.__fc.setTolerance(tolerance)
337
338 def __reinitialise(self,phi):
339 print "reintialization started:"
340 s=self.__makeInterface(phi,1.)
341 print "phi range:",inf(phi), sup(phi)
342 g=grad(phi)
343 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 # while (iter<=self.__reinit_max):
354 # 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 # while (iter<=self.__reinit_max):
367 # 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 if self.__reinitFC:
377 self.__reinitfc.setValue(C=-w,Y=s)
378 self.__reinitfc.setInitialSolution(phi+self.__diam)
379 dtau=self.__reinitfc.getSafeTimeStepSize()
380 # dtau = 0.5*inf(Function(self.__domain).getSize())
381 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 iter =0
385 # self.__reinitPde.setValue(q=whereNegative(abs(phi)-2*self.__h), r=phi)
386 # self.__reinitPde.setValue(r=phi)
387 while (iter<=self.__reinit_max):
388 phi_old=phi
389 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 change = Lsup(phi-phi_old)/self.__diam
400 print "phi range:",inf(phi), sup(phi)
401 print "iteration :", iter, " change:", change
402 iter +=1
403 print "reintialization completed."
404 return phi
405
406 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
410
411 #================ 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 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 def __reinitialise_old(self):
440 #=============================================
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 while c < self.__reinit_max: # and abs(diff) >= 0.01:
470 #
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 # self.__reinitPde.setValue(Y =self.__phi, Y_reduced=dtau*s*(1.-inner(w,grad(phi_half,ReducedFunction(self.__domain)))))
482 # 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 # self.__updateInterface()
489
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 x2=self.__domain.getX()
552 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