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

  ViewVC Help
Powered by ViewVC 1.1.26