/[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 1716 - (hide annotations)
Thu Aug 21 05:03:49 2008 UTC (11 years, 6 months ago) by gross
File MIME type: text/x-python
File size: 14966 byte(s)
getListOfTags method added to FunctionSpace class
1 gross 1473 from esys.escript import *
2 gross 1659 from esys.escript.linearPDEs import LinearPDE, TransportPDE
3 gross 1473 from esys.escript.pdetools import Projector
4     from esys import finley
5     import math
6    
7    
8    
9     class LevelSet(object):
10 gross 1659 def __init__(self,phi,reinitialization_steps_max=10,reinitialization_each=2,relative_smoothing_width=2.):
11 gross 1473 """
12     initialize model
13     """
14 gross 1659 self.__domain = phi.getDomain()
15     x=self.__domain.getX()
16     diam=0
17     for i in range(self.__domain.getDim()):
18     xi=x[i]
19     diam+=(inf(xi)-sup(xi))**2
20     self.__diam=sqrt(diam)
21 gross 1661 self.__h = sup(Function(self.__domain).getSize())
22 gross 1659 self.__reinitialization_each=reinitialization_each
23     self.__reinitialization_steps_max = reinitialization_steps_max
24     self.__relative_smoothing_width = relative_smoothing_width
25     self.__phi = phi
26     self.__update_count=0
27     self.velocity = None
28 gross 1473
29 gross 1659 #==================================================
30 gross 1673 self.__FC=True
31     if self.__FC:
32     self.__fc=TransportPDE(self.__domain,num_equations=1,theta=0.5)
33     # self.__fc.setReducedOrderOn()
34     self.__fc.setValue(M=Scalar(1.,Function(self.__domain)))
35     self.__fc.setInitialSolution(phi+self.__diam)
36     else:
37     self.__fcpde = LinearPDE(self.__domain,numEquations=1, numSolutions=1)
38     self.__fcpde.setSolverMethod(solver=LinearPDE.LUMPING)
39     self.__fcpde.setReducedOrderOn()
40     self.__fcpde.setSymmetryOn()
41     self.__fcpde.setValue(D=1.)
42 gross 1659
43     #=======================================
44 gross 1673 self.__reinitFC=False
45     if self.__reinitFC:
46 gross 1716 self.__reinitfc=TransportPDE(self.__domain,num_equations=1,theta=1.0)
47 gross 1673 self.__reinitfc.setValue(M=Scalar(1.,Function(self.__domain)))
48     self.__reinitfc.setTolerance(1.e-5)
49     else:
50     self.__reinitPde = LinearPDE(self.__domain,numEquations=1, numSolutions=1)
51     # self.__reinitPde.setSolverMethod(solver=LinearPDE.LUMPING)
52     self.__reinitPde.setReducedOrderOn()
53     self.__reinitPde.setSymmetryOn()
54     self.__reinitPde.setValue(D=1.)
55 gross 1716 self.__reinitPde.setTolerance(1.e-2)
56 gross 1673 # self.__reinitPde.setSolverMethod(preconditioner=LinearPDE.ILU0)
57 gross 1659 #=======================================
58 gross 1661 self.__updateInterface()
59     print "phi range:",inf(phi), sup(phi)
60 gross 1473
61 gross 1661 def setFixedRegion(self,mask,contour=0):
62     q=whereNonPositive(abs(self.__phi-contour)-self.__h)*mask
63     # self.__fc.setValue(q=q)
64     self.__reinitPde.setValue(q=q)
65 gross 1473
66 gross 1673 def getH(self):
67     return self.__h
68 gross 1473 def getDomain(self):
69 gross 1659 return self.__domain
70 gross 1473
71 gross 1659 def getTimeStepSize(self,velocity):
72     """
73     returns a new dt for a given velocity using the courant coundition
74     """
75     self.velocity=velocity
76 gross 1673 if self.__FC:
77     self.__fc.setValue(C=-interpolate(velocity,Function(self.__domain)))
78     dt=self.__fc.getSafeTimeStepSize()
79     else:
80     dt=0.5*self.__h/sup(length(velocity))
81 gross 1659 return dt
82    
83     def getLevelSetFunction(self):
84     """
85     returns the level set function
86     """
87     return self.__phi
88    
89 gross 1661 def __updateInterface(self):
90     self.__smoothed_char=self.__makeInterface(self.__phi)
91    
92     def __makeInterface(self,phi,smoothing_width=1.):
93 gross 1659 """
94     creates a very smooth interface from -1 to 1 over the length 2*h*smoothing_width where -1 is used where the level set is negative
95     and 1 where the level set is 1
96     """
97 gross 1473 s=smoothing_width*self.__h
98 gross 1661 phi_on_h=interpolate(phi,Function(self.__domain))
99 gross 1473 mask_neg = whereNonNegative(-s-phi_on_h)
100     mask_pos = whereNonNegative(phi_on_h-s)
101     mask_interface = 1.-mask_neg-mask_pos
102 gross 1661 # interface=1.-(phi_on_h-s)**2/(2.*s**3)*(phi_on_h+2*s) # function f with f(s)=1, f'(s)=0, f(-s)=-1, f'(-s)=0, f(0)=0, f'(0)=
103     interface=phi_on_h/s
104 gross 1473 return - mask_neg + mask_pos + mask_interface * interface
105    
106 gross 1661 def makeCharacteristicFunction(self, contour=0, positiveSide=True, smoothing_width=1.):
107     return self.makeCharacteristicFunctionFromExternalLevelSetFunction(self.__phi,contour,positiveSide,smoothing_width)
108    
109     def makeCharacteristicFunctionFromExternalLevelSetFunction(self,phi,contour=0, positiveSide=True, smoothing_width=1.):
110     s=self.__makeInterface(phi-contour,smoothing_width)
111     if positiveSide:
112     return (1+s)/2
113     else:
114     return (1-s)/2
115    
116    
117 gross 1473 def update(self,dt):
118     """
119     sets a new velocity and updates the level set fuction
120    
121     @param dt: time step forward
122     """
123 gross 1673 if self.__FC:
124     self.__phi=self.__fc.solve(dt, verbose=False)-self.__diam
125     else:
126     self.__fcpde.setValue(Y = self.__phi-(dt/2.)*inner(self.velocity,grad(self.__phi)))
127     phi_half = self.__fcpde.getSolution()
128     self.__fcpde.setValue(Y = self.__phi-dt*inner(self.velocity,grad(phi_half)))
129     self.__phi= self.__fcpde.getSolution()
130 gross 1659 self.__update_count += 1
131 gross 1473 if self.__update_count%self.__reinitialization_each == 0:
132 gross 1661 self.__phi=self.__reinitialise(self.__phi)
133 gross 1673 if self.__FC:
134     self.__fc.setInitialSolution(self.__phi+self.__diam)
135 gross 1661 self.__updateInterface()
136 gross 1473
137 gross 1659 def setTolerance(self,tolerance=1e-3):
138 gross 1673 if self.__FC:
139     self.__fc.setTolerance(tolerance)
140 gross 1473
141 gross 1659 def __reinitialise(self,phi):
142     print "reintialization started:"
143 gross 1716 s=self.__makeInterface(phi,1.)
144     print "phi range:",inf(phi), sup(phi)
145 gross 1659 g=grad(phi)
146 gross 1716 w=s*g/(length(g)+EPSILON)
147     #=======================================================
148     # # positive part:
149     # phi_p=wherePositive(phi)*phi
150     # self.__reinitfc.setInitialSolution(phi_p)
151     # self.__reinitfc.setValue(C=-wherePositive(s)*w,Y=wherePositive(s)*s,q=whereNonPositive(phi))
152     # dtau=self.__h
153     # print "step size: dt (pos)= ",dtau
154     # print "phi_p range:",inf(phi_p), sup(phi_p)
155     # iter=0
156     # while (iter<=self.__reinitialization_steps_max):
157     # phi_p=self.__reinitfc.solve(dtau)
158     # print "phi_p range:",inf(phi_p), sup(phi_p)
159     # iter+=1
160     # # negative part:
161     # phi_n=-whereNegative(phi)*phi
162     # self.__reinitfc.setInitialSolution(phi_n)
163     # self.__reinitfc.setValue(C=-whereNegative(s)*w,Y=-whereNegative(s)*s,q=whereNonNegative(phi))
164     # # dtau=self.__reinitfc.getSafeTimeStepSize()
165     # dtau=self.__h
166     # print "step size: dt (neg)= ",dtau
167     # print "phi_n range:",inf(phi_n), sup(phi_n)
168     # iter=0
169     # while (iter<=self.__reinitialization_steps_max):
170     # phi_n=self.__reinitfc.solve(dtau)
171     # print "phi_n range:",inf(phi_n), sup(phi_n)
172     # iter+=1
173     # phi=phi_p-phi_n
174     # print "phi range:",inf(phi), sup(phi)
175     # print "reintialization completed."
176     # return phi
177    
178     #=======================================================
179 gross 1673 if self.__reinitFC:
180     self.__reinitfc.setValue(C=-w,Y=s)
181     self.__reinitfc.setInitialSolution(phi+self.__diam)
182 gross 1716 dtau=self.__reinitfc.getSafeTimeStepSize()
183     # dtau = 0.5*inf(Function(self.__domain).getSize())
184 gross 1673 else:
185     dtau = 0.5*inf(Function(self.__domain).getSize())
186     print "step size: dt = ",dtau,inf(abs(phi.interpolate(Function(self.__domain)))/abs(s-inner(w,grad(phi))))
187 gross 1659 iter =0
188 gross 1673 # self.__reinitPde.setValue(q=whereNegative(abs(phi)-2*self.__h), r=phi)
189 gross 1661 # self.__reinitPde.setValue(r=phi)
190 gross 1659 while (iter<=self.__reinitialization_steps_max):
191     phi_old=phi
192 gross 1673 if self.__reinitFC:
193     phi = self.__reinitfc.solve(dtau)-self.__diam
194     else:
195     self.__reinitPde.setValue(Y = phi+(dtau/2.)*(s-inner(w,grad(phi))))
196     phi_half = self.__reinitPde.getSolution()
197     self.__reinitPde.setValue(Y = phi+dtau*(s-inner(w,grad(phi_half))))
198     # g=grad(phi)
199     # S=inner(w,grad(phi))
200     # self.__reinitPde.setValue(Y = phi+dtau*(s-S),X=dtau**2/2*w*(s-S))
201     phi = self.__reinitPde.getSolution()
202 gross 1659 change = Lsup(phi-phi_old)/self.__diam
203     print "phi range:",inf(phi), sup(phi)
204 gross 1716 print "iteration :", iter, " change:", change
205 gross 1659 iter +=1
206     print "reintialization completed."
207     return phi
208 gross 1473
209 gross 1661 def createParameter(self,value_negative=-1.,value_positive=1):
210     out = (value_negative+value_positive)/2. + self.__smoothed_char * ((value_positive-value_negative)/2.)
211     return out
212 gross 1473
213 gross 1661
214 gross 1659 #================ things from here onwards are not used nor tested: ==========================================
215    
216    
217     def getCharacteristicFunction(self):
218     return self.__smoothed_char
219    
220    
221    
222    
223 gross 1473 def RK2(self,L):
224     k0=L(phi)
225     phi_1=phi+dt*k0
226     k1=L(phi_1)
227     phi=phi+dt/2*(k0+k1)
228     def RK3(self,L):
229     k0=L(phi)
230     phi_1=phi+dt*L(phi)
231     k1=L(phi_1)
232     phi_2=phi+dt/4*(k0+k1)
233     k2=L(phi_2)
234     phi=phi+dt/6*(k0+4*k2+K1)
235     def TG(self,L):
236     k0=L(phi)
237     phi_1=phi+dt/2*k0
238     k1=L(phi_1)
239     phi=phi+dt*k1
240    
241    
242 gross 1659 def __reinitialise_old(self):
243 gross 1473 #=============================================
244     f=0.1
245     f2=0.3
246     h=self.__h
247     fs=h.getFunctionSpace()
248     vol=self.getVolumeOfNegativeDomain()
249     self.__reinitPde.setValue(D=1.0)
250     #============
251     grad_phi=grad(self.__phi,fs)
252     len_grad_phi=length(grad_phi)
253     w=grad_phi/len_grad_phi
254    
255     # s=wherePositive(self. __makeInterface(2.))-whereNegative(self. __makeInterface(2.))
256     s=self. __makeInterface(2.)
257     # s2=whereNegative(len_grad_phi-f2)*(1-(1-len_grad_phi/f2)**2)+(1-whereNegative(len_grad_phi-f2))
258     s2=1.
259     s=s*s2
260     dtau=f*inf(h/abs(s))
261    
262     #========================
263     diff=1.
264     d=1.
265     c =0
266     #==============================================
267     TVD=integrate(length(grad_phi))
268     print "initial range ",inf(self.__phi),sup(self.__phi),"error:",Lsup(1.-len_grad_phi),"volume =",vol,TVD
269     # saveVTK("test.%s.xml"%c,l=length(grad(self.__phi,fs))-1,s=s,phi=self.__phi)
270    
271     dtau=f*inf(h/abs(s))
272     while c < self.__reinitialization_steps_max: # and abs(diff) >= 0.01:
273     #
274     grad_phi=grad(self.__phi,fs)
275     len_grad_phi=length(grad_phi)
276     #
277     # s=self.__makeInterface(1.)
278     # s2=whereNegative(len_grad_phi-f2)*(1-(1-len_grad_phi/f2)**2)+(1-whereNegative(len_grad_phi-f2))
279     # s=s*s2
280     # phi_on_h=interpolate(self.__phi,fs)
281     # self.__reinitPde.setValue(Y =self.__phi+dtau/2*s*(1.-len_grad_phi))
282     # phi_half=self.__reinitPde.getSolution()
283     # self.__reinitPde.setValue(Y =self.__phi+dtau*s*(1.-length(grad(phi_half,fs))))
284 gross 1659 # self.__reinitPde.setValue(Y =self.__phi, Y_reduced=dtau*s*(1.-inner(w,grad(phi_half,ReducedFunction(self.__domain)))))
285 gross 1473 # self.__reinitPde.setValue(Y =self.__phi+dtau*(s-inner(w,grad_phi)),X=dtau/2*h*(s-inner(w,grad_phi))*w)
286     # self.__reinitPde.setValue(Y =self.__phi+dtau*s*(1.-len_grad_phi),X=-dtau/2*h*s**2*grad_phi)
287     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)
288    
289     # 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)
290     self.__phi, previous = self.__reinitPde.getSolution(), self.__phi
291 gross 1661 # self.__updateInterface()
292 gross 1473
293     vol,vol_old=self.getVolumeOfNegativeDomain(),vol
294     diff=(vol-vol_old)
295     r=Lsup(length(grad(self.__phi))-1.)
296     TVD=integrate(length(grad(self.__phi,fs)))
297     print "iteration :", c, "range ",inf(self.__phi),sup(self.__phi),"error :",r,"volume change:",diff,TVD
298     # saveVTK("test.%s.xml"%(c+1),l=length(grad(self.__phi,fs)),s=s,phi=self.__phi,v=grad(self.__phi,fs))
299     c += 1
300     return
301     #==============================================
302     f2=0.7
303     h=self.__h
304     fs=h.getFunctionSpace()
305     vol=self.getVolumeOfNegativeDomain()
306     s=abs(self. __makeInterface(2.))
307     grad_phi=grad(self.__phi,fs)
308     len_grad_phi=length(grad_phi)
309    
310     #----
311     # aphi_on_h=abs(interpolate(self.__phi,self.__h.getFunctionSpace()))
312     # q=2*self.__h
313     # mask_neg = whereNegative(aphi_on_h-q)
314     # mask_pos = wherePositive(aphi_on_h-q-self.__h)
315     # mask_interface = 1.-mask_neg-mask_pos
316     # s2=mask_neg + mask_interface * (q+self.__h-aphi_on_h)/self.__h
317     #----
318     m=whereNonPositive(len_grad_phi-f2)
319     s2=m*len_grad_phi/f2+(1.-m)
320     #----
321     # s2=1.
322     #----
323     self.__reinitPde.setValue(D=(1-s)/h**2,Y=(1-s)/h**2*self.__phi)
324     self.__reinitPde.setValue(A=(s2*len_grad_phi+(1.-s2))*kronecker(3),X=grad_phi)
325     # self.__reinitPde.setValue(A=(len_grad_phi-1)/len_grad_phi*kronecker(3))
326     # self.__reinitPde.setValue(A=kronecker(3), X=grad_phi/len_grad_phi)
327     self.__phi = self.__reinitPde.getSolution()
328     r=Lsup(length(grad(self.__phi))-1.)
329     vol,vol_old=self.getVolumeOfNegativeDomain(),vol
330     diff=(vol-vol_old)/vol
331     print "iteration :", inf(self.__phi),sup(self.__phi),r,diff
332     # saveVTK("test.%s.xml"%0,l=length(grad(self.__phi,fs)),s=s,phi=self.__phi,v=grad(self.__phi,fs),s2=s2)
333     return
334     #=============================================
335    
336     #============
337    
338    
339     def getVolumeOfNegativeDomain(self):
340     """
341     return the current volume of domain with phi<0.
342     """
343     return integrate((1.-self.__makeInterface(1.))/2.)
344    
345    
346     def getBoundingBoxOfNegativeDomain(self):
347     """
348     get the height of the region with phi<0
349     """
350     fs=self.__h.getFunctionSpace()
351     mask_phi1=wherePositive(interpolate(self.__phi,fs))
352     mask_phi2=wherePositive(self.__phi)
353     x1=fs.getX()
354 gross 1659 x2=self.__domain.getX()
355 gross 1473 out=[]
356     for i in range(fs.getDim()):
357     x1_i=x1[i]
358     x2_i=x2[i]
359     d=2*(sup(x2_i)-inf(x2_i))
360     offset1=d*mask_phi1
361     offset2=d*mask_phi2
362     out.append((min(inf(x1_i+offset1),inf(x2_i+offset2)),max(sup(x1_i-offset1),sup(x2_i-offset2))))
363     return tuple(out)
364    

  ViewVC Help
Powered by ViewVC 1.1.26