/[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 2158 - (hide annotations)
Mon Dec 15 07:17:47 2008 UTC (11 years, 2 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 ksteube 1811
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 gross 1473 from esys.escript import *
23 gross 1659 from esys.escript.linearPDEs import LinearPDE, TransportPDE
24 gross 1473 from esys.escript.pdetools import Projector
25     from esys import finley
26     import math
27    
28    
29    
30 gross 1788 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 caltinay 2158 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 gross 1788
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 caltinay 2158
70 gross 1788 def __advect(self, velocity, dt):
71     """
72 caltinay 2158 advects the level set function in the presence of a velocity field.
73 gross 1788
74 caltinay 2158 This implementation uses the 2-step Taylor-Galerkin method
75 gross 1788 @param velocity: velocity field
76 caltinay 2158 @param dt: time increment
77 gross 1788 @return: the advected level set function
78     """
79     Y = self.__phi-dt/2.*inner(velocity,grad(self.__phi))
80 caltinay 2158 self.__PDE.setValue(Y=Y)
81 gross 1788 phi_half = self.__PDE.getSolution()
82     Y = self.__phi-dt*inner(velocity,grad(phi_half))
83 caltinay 2158 self.__PDE.setValue(Y=Y)
84 gross 1788 phi = self.__PDE.getSolution()
85     print "LevelSet: Advection step done"
86     return phi
87    
88     def __reinitialise(self):
89     """
90 caltinay 2158 reinitializes the level set
91 gross 1788
92 caltinay 2158 It solves the...
93 gross 1788
94 caltinay 2158 @return: reinitialized level set
95 gross 1788 """
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 caltinay 2158 returns a new dt for a given velocity using the courant condition
121 gross 1788
122     @param velocity: velocity field
123     """
124     self.__velocity=velocity
125     dt=0.5*self.__h/sup(length(velocity))
126     return dt
127 caltinay 2158
128 gross 1788 def update(self,dt):
129     """
130 caltinay 2158 sets a new velocity and updates the level set function
131 gross 1788
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 caltinay 2158 updates phi under the presence of a velocity field
142 gross 1788
143 caltinay 2158 If dt is small this call is equivalent to call
144    
145 gross 1788 dt=LevelSet.getTimeStepSize(velocity)
146     phi=LevelSet.update(dt)
147 caltinay 2158
148 gross 1788 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 caltinay 2158 returns the volume of the phi(x)<0 region
164 gross 1788 """
165     return integrate(whereNegative(self.__phi.interpolate(Function(self.__domain))))
166    
167     def getSurface(self,rel_width_factor=0.5):
168     """
169 caltinay 2158 returns a mask for the phi(x)=1 region
170    
171     @param rel_width_factor: relative width of region around zero contour.
172 gross 1788 """
173     return whereNegative(abs(self.__phi)-rel_width_factor*self.__h)
174 caltinay 2158
175 gross 1788 def getH(self):
176     """
177 caltinay 2158 returns the mesh size
178 gross 1788 """
179     return self.__h
180    
181     def getDomain(self):
182     """
183 caltinay 2158 returns the domain
184 gross 1788 """
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 caltinay 2158 creates a function with param_neg where phi<0 and param_pos where phi>0
196     (no smoothing)
197 gross 1788
198     @param param_neg: value of parameter on the negative side (phi<0)
199 caltinay 2158 @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 gross 1788 """
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 caltinay 2158 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 gross 1788
212 caltinay 2158 @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 gross 1788 """
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 caltinay 2158 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 gross 1788 """
226     s=smoothing_width*self.__h
227 caltinay 2158 phi_on_h=interpolate(phi,Function(self.__domain))
228 gross 1788 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 caltinay 2158 makes a smooth characteristic function of the region phi(x)>contour if
237     positiveSide and phi(x)<contour otherwise.
238 gross 1788
239 caltinay 2158 @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 gross 1788 """
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 gross 1473 """
259     initialize model
260     """
261 caltinay 2158 self.__domain = phi.getDomain()
262 gross 1659 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 gross 1661 self.__h = sup(Function(self.__domain).getSize())
269 gross 1788 self.__reinit_each=reinit_each
270     self.__reinit_max = reinit_max
271     self.__smooth = smooth
272 gross 1659 self.__phi = phi
273     self.__update_count=0
274     self.velocity = None
275 gross 1473
276 gross 1659 #==================================================
277 gross 1673 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 caltinay 2158 self.__fcpde.setSolverMethod(solver=LinearPDE.LUMPING)
286 gross 1673 self.__fcpde.setReducedOrderOn()
287     self.__fcpde.setSymmetryOn()
288     self.__fcpde.setValue(D=1.)
289 gross 1659
290     #=======================================
291 gross 1673 self.__reinitFC=False
292     if self.__reinitFC:
293 gross 1716 self.__reinitfc=TransportPDE(self.__domain,num_equations=1,theta=1.0)
294 gross 1673 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 caltinay 2158 # self.__reinitPde.setSolverMethod(solver=LinearPDE.LUMPING)
299 gross 1673 self.__reinitPde.setReducedOrderOn()
300     self.__reinitPde.setSymmetryOn()
301     self.__reinitPde.setValue(D=1.)
302 gross 1716 self.__reinitPde.setTolerance(1.e-2)
303 gross 1673 # self.__reinitPde.setSolverMethod(preconditioner=LinearPDE.ILU0)
304 gross 1659 #=======================================
305 gross 1661 self.__updateInterface()
306     print "phi range:",inf(phi), sup(phi)
307 gross 1473
308 gross 1661 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 gross 1473
313    
314 gross 1661 def __updateInterface(self):
315     self.__smoothed_char=self.__makeInterface(self.__phi)
316    
317 gross 1473
318     def update(self,dt):
319     """
320 caltinay 2158 sets a new velocity and updates the level set function
321 gross 1473
322     @param dt: time step forward
323     """
324 gross 1673 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 gross 1659 self.__update_count += 1
332 gross 1788 if self.__update_count%self.__reinit_each == 0:
333 gross 1661 self.__phi=self.__reinitialise(self.__phi)
334 gross 1673 if self.__FC:
335     self.__fc.setInitialSolution(self.__phi+self.__diam)
336 gross 1661 self.__updateInterface()
337 gross 1473
338 gross 1659 def setTolerance(self,tolerance=1e-3):
339 gross 1673 if self.__FC:
340     self.__fc.setTolerance(tolerance)
341 gross 1473
342 gross 1659 def __reinitialise(self,phi):
343     print "reintialization started:"
344 gross 1716 s=self.__makeInterface(phi,1.)
345     print "phi range:",inf(phi), sup(phi)
346 gross 1659 g=grad(phi)
347 gross 1716 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 gross 1788 # while (iter<=self.__reinit_max):
358 gross 1716 # 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 caltinay 2158 # # dtau=self.__reinitfc.getSafeTimeStepSize()
366 gross 1716 # 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 gross 1788 # while (iter<=self.__reinit_max):
371 gross 1716 # 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 gross 1673 if self.__reinitFC:
381     self.__reinitfc.setValue(C=-w,Y=s)
382     self.__reinitfc.setInitialSolution(phi+self.__diam)
383 caltinay 2158 dtau=self.__reinitfc.getSafeTimeStepSize()
384 gross 1716 # dtau = 0.5*inf(Function(self.__domain).getSize())
385 gross 1673 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 gross 1659 iter =0
389 gross 1673 # self.__reinitPde.setValue(q=whereNegative(abs(phi)-2*self.__h), r=phi)
390 gross 1661 # self.__reinitPde.setValue(r=phi)
391 gross 1788 while (iter<=self.__reinit_max):
392 gross 1659 phi_old=phi
393 gross 1673 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 gross 1659 change = Lsup(phi-phi_old)/self.__diam
404     print "phi range:",inf(phi), sup(phi)
405 gross 1716 print "iteration :", iter, " change:", change
406 gross 1659 iter +=1
407     print "reintialization completed."
408     return phi
409 gross 1473
410 gross 1661 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 gross 1473
414 gross 1661
415 caltinay 2158 #=========== things from here onwards are not used nor tested: ===========
416 gross 1659
417 caltinay 2158
418 gross 1659 def getCharacteristicFunction(self):
419     return self.__smoothed_char
420    
421    
422 gross 1473 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 caltinay 2158
428 gross 1473 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 caltinay 2158
436 gross 1473 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 caltinay 2158
442    
443 gross 1659 def __reinitialise_old(self):
444 gross 1473 #=============================================
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 gross 1788 while c < self.__reinit_max: # and abs(diff) >= 0.01:
474 gross 1473 #
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 caltinay 2158 # phi_on_h=interpolate(self.__phi,fs)
482 gross 1473 # 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 gross 1659 # self.__reinitPde.setValue(Y =self.__phi, Y_reduced=dtau*s*(1.-inner(w,grad(phi_half,ReducedFunction(self.__domain)))))
486 gross 1473 # 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 gross 1661 # self.__updateInterface()
493 gross 1473
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 caltinay 2158
511 gross 1473 #----
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 caltinay 2158 returns the current volume of domain with phi<0.
543 gross 1473 """
544     return integrate((1.-self.__makeInterface(1.))/2.)
545    
546 caltinay 2158
547 gross 1473 def getBoundingBoxOfNegativeDomain(self):
548     """
549 caltinay 2158 returns the height of the region with phi<0
550 gross 1473 """
551     fs=self.__h.getFunctionSpace()
552     mask_phi1=wherePositive(interpolate(self.__phi,fs))
553     mask_phi2=wherePositive(self.__phi)
554     x1=fs.getX()
555 gross 1659 x2=self.__domain.getX()
556 gross 1473 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