/[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 4114 - (show annotations)
Fri Dec 14 04:24:46 2012 UTC (6 years, 8 months ago) by caltinay
File MIME type: text/x-python
File size: 21506 byte(s)
Time to remove deprecated saveVTK/DX methods from Data and Domain.

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

  ViewVC Help
Powered by ViewVC 1.1.26