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

  ViewVC Help
Powered by ViewVC 1.1.26