/[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 2158 - (show annotations)
Mon Dec 15 07:17:47 2008 UTC (10 years, 10 months ago) by caltinay
File MIME type: text/x-python
File size: 20862 byte(s)
Assorted spelling, grammar, whitespace and copy/paste error fixes (Part 1).
This commit should be a no-op.

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

  ViewVC Help
Powered by ViewVC 1.1.26