/[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 1673 - (show annotations)
Thu Jul 24 22:28:50 2008 UTC (11 years, 2 months ago) by gross
File MIME type: text/x-python
File size: 13490 byte(s)
more algorithms for level set
1 from esys.escript import *
2 from esys.escript.linearPDEs import LinearPDE, TransportPDE
3 from esys.escript.pdetools import Projector
4 from esys import finley
5 import math
6
7
8
9 class LevelSet(object):
10 def __init__(self,phi,reinitialization_steps_max=10,reinitialization_each=2,relative_smoothing_width=2.):
11 """
12 initialize model
13 """
14 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 self.__h = sup(Function(self.__domain).getSize())
22 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
29 #==================================================
30 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
43 #=======================================
44 self.__reinitFC=False
45 if self.__reinitFC:
46 self.__reinitfc=TransportPDE(self.__domain,num_equations=1,theta=1.)
47 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 # self.__reinitPde.setTolerance(1.e-8)
56 # self.__reinitPde.setSolverMethod(preconditioner=LinearPDE.ILU0)
57 #=======================================
58 self.__updateInterface()
59 print "phi range:",inf(phi), sup(phi)
60
61 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
66 def getH(self):
67 return self.__h
68 def getDomain(self):
69 return self.__domain
70
71 def getTimeStepSize(self,velocity):
72 """
73 returns a new dt for a given velocity using the courant coundition
74 """
75 self.velocity=velocity
76 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 return dt
82
83 def getLevelSetFunction(self):
84 """
85 returns the level set function
86 """
87 return self.__phi
88
89 def __updateInterface(self):
90 self.__smoothed_char=self.__makeInterface(self.__phi)
91
92 def __makeInterface(self,phi,smoothing_width=1.):
93 """
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 s=smoothing_width*self.__h
98 phi_on_h=interpolate(phi,Function(self.__domain))
99 mask_neg = whereNonNegative(-s-phi_on_h)
100 mask_pos = whereNonNegative(phi_on_h-s)
101 mask_interface = 1.-mask_neg-mask_pos
102 # 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 return - mask_neg + mask_pos + mask_interface * interface
105
106 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 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 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 self.__update_count += 1
131 if self.__update_count%self.__reinitialization_each == 0:
132 self.__phi=self.__reinitialise(self.__phi)
133 if self.__FC:
134 self.__fc.setInitialSolution(self.__phi+self.__diam)
135 self.__updateInterface()
136
137 def setTolerance(self,tolerance=1e-3):
138 if self.__FC:
139 self.__fc.setTolerance(tolerance)
140
141 def __reinitialise(self,phi):
142 print "reintialization started:"
143 s=self.__makeInterface(phi,1.5)
144 g=grad(phi)
145 w = s*g/(length(g)+EPSILON)
146 if self.__reinitFC:
147 self.__reinitfc.setValue(C=-w,Y=s)
148 self.__reinitfc.setInitialSolution(phi+self.__diam)
149 # dtau=self.__reinitfc.getSafeTimeStepSize()
150 dtau = 0.5*inf(Function(self.__domain).getSize())
151 else:
152 dtau = 0.5*inf(Function(self.__domain).getSize())
153 print "step size: dt = ",dtau,inf(abs(phi.interpolate(Function(self.__domain)))/abs(s-inner(w,grad(phi))))
154 iter =0
155 # self.__reinitPde.setValue(q=whereNegative(abs(phi)-2*self.__h), r=phi)
156 # self.__reinitPde.setValue(r=phi)
157 while (iter<=self.__reinitialization_steps_max):
158 phi_old=phi
159 if self.__reinitFC:
160 phi = self.__reinitfc.solve(dtau)-self.__diam
161 else:
162 self.__reinitPde.setValue(Y = phi+(dtau/2.)*(s-inner(w,grad(phi))))
163 phi_half = self.__reinitPde.getSolution()
164 self.__reinitPde.setValue(Y = phi+dtau*(s-inner(w,grad(phi_half))))
165 # g=grad(phi)
166 # S=inner(w,grad(phi))
167 # self.__reinitPde.setValue(Y = phi+dtau*(s-S),X=dtau**2/2*w*(s-S))
168 phi = self.__reinitPde.getSolution()
169 change = Lsup(phi-phi_old)/self.__diam
170 print "phi range:",inf(phi), sup(phi)
171 print "iteration :", iter, " error:", change
172 iter +=1
173 print "reintialization completed."
174 return phi
175
176 def createParameter(self,value_negative=-1.,value_positive=1):
177 out = (value_negative+value_positive)/2. + self.__smoothed_char * ((value_positive-value_negative)/2.)
178 return out
179
180
181 #================ things from here onwards are not used nor tested: ==========================================
182
183
184 def getCharacteristicFunction(self):
185 return self.__smoothed_char
186
187
188
189
190 def RK2(self,L):
191 k0=L(phi)
192 phi_1=phi+dt*k0
193 k1=L(phi_1)
194 phi=phi+dt/2*(k0+k1)
195 def RK3(self,L):
196 k0=L(phi)
197 phi_1=phi+dt*L(phi)
198 k1=L(phi_1)
199 phi_2=phi+dt/4*(k0+k1)
200 k2=L(phi_2)
201 phi=phi+dt/6*(k0+4*k2+K1)
202 def TG(self,L):
203 k0=L(phi)
204 phi_1=phi+dt/2*k0
205 k1=L(phi_1)
206 phi=phi+dt*k1
207
208
209 def __reinitialise_old(self):
210 #=============================================
211 f=0.1
212 f2=0.3
213 h=self.__h
214 fs=h.getFunctionSpace()
215 vol=self.getVolumeOfNegativeDomain()
216 self.__reinitPde.setValue(D=1.0)
217 #============
218 grad_phi=grad(self.__phi,fs)
219 len_grad_phi=length(grad_phi)
220 w=grad_phi/len_grad_phi
221
222 # s=wherePositive(self. __makeInterface(2.))-whereNegative(self. __makeInterface(2.))
223 s=self. __makeInterface(2.)
224 # s2=whereNegative(len_grad_phi-f2)*(1-(1-len_grad_phi/f2)**2)+(1-whereNegative(len_grad_phi-f2))
225 s2=1.
226 s=s*s2
227 dtau=f*inf(h/abs(s))
228
229 #========================
230 diff=1.
231 d=1.
232 c =0
233 #==============================================
234 TVD=integrate(length(grad_phi))
235 print "initial range ",inf(self.__phi),sup(self.__phi),"error:",Lsup(1.-len_grad_phi),"volume =",vol,TVD
236 # saveVTK("test.%s.xml"%c,l=length(grad(self.__phi,fs))-1,s=s,phi=self.__phi)
237
238 dtau=f*inf(h/abs(s))
239 while c < self.__reinitialization_steps_max: # and abs(diff) >= 0.01:
240 #
241 grad_phi=grad(self.__phi,fs)
242 len_grad_phi=length(grad_phi)
243 #
244 # s=self.__makeInterface(1.)
245 # s2=whereNegative(len_grad_phi-f2)*(1-(1-len_grad_phi/f2)**2)+(1-whereNegative(len_grad_phi-f2))
246 # s=s*s2
247 # phi_on_h=interpolate(self.__phi,fs)
248 # self.__reinitPde.setValue(Y =self.__phi+dtau/2*s*(1.-len_grad_phi))
249 # phi_half=self.__reinitPde.getSolution()
250 # self.__reinitPde.setValue(Y =self.__phi+dtau*s*(1.-length(grad(phi_half,fs))))
251 # self.__reinitPde.setValue(Y =self.__phi, Y_reduced=dtau*s*(1.-inner(w,grad(phi_half,ReducedFunction(self.__domain)))))
252 # self.__reinitPde.setValue(Y =self.__phi+dtau*(s-inner(w,grad_phi)),X=dtau/2*h*(s-inner(w,grad_phi))*w)
253 # self.__reinitPde.setValue(Y =self.__phi+dtau*s*(1.-len_grad_phi),X=-dtau/2*h*s**2*grad_phi)
254 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)
255
256 # 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)
257 self.__phi, previous = self.__reinitPde.getSolution(), self.__phi
258 # self.__updateInterface()
259
260 vol,vol_old=self.getVolumeOfNegativeDomain(),vol
261 diff=(vol-vol_old)
262 r=Lsup(length(grad(self.__phi))-1.)
263 TVD=integrate(length(grad(self.__phi,fs)))
264 print "iteration :", c, "range ",inf(self.__phi),sup(self.__phi),"error :",r,"volume change:",diff,TVD
265 # saveVTK("test.%s.xml"%(c+1),l=length(grad(self.__phi,fs)),s=s,phi=self.__phi,v=grad(self.__phi,fs))
266 c += 1
267 return
268 #==============================================
269 f2=0.7
270 h=self.__h
271 fs=h.getFunctionSpace()
272 vol=self.getVolumeOfNegativeDomain()
273 s=abs(self. __makeInterface(2.))
274 grad_phi=grad(self.__phi,fs)
275 len_grad_phi=length(grad_phi)
276
277 #----
278 # aphi_on_h=abs(interpolate(self.__phi,self.__h.getFunctionSpace()))
279 # q=2*self.__h
280 # mask_neg = whereNegative(aphi_on_h-q)
281 # mask_pos = wherePositive(aphi_on_h-q-self.__h)
282 # mask_interface = 1.-mask_neg-mask_pos
283 # s2=mask_neg + mask_interface * (q+self.__h-aphi_on_h)/self.__h
284 #----
285 m=whereNonPositive(len_grad_phi-f2)
286 s2=m*len_grad_phi/f2+(1.-m)
287 #----
288 # s2=1.
289 #----
290 self.__reinitPde.setValue(D=(1-s)/h**2,Y=(1-s)/h**2*self.__phi)
291 self.__reinitPde.setValue(A=(s2*len_grad_phi+(1.-s2))*kronecker(3),X=grad_phi)
292 # self.__reinitPde.setValue(A=(len_grad_phi-1)/len_grad_phi*kronecker(3))
293 # self.__reinitPde.setValue(A=kronecker(3), X=grad_phi/len_grad_phi)
294 self.__phi = self.__reinitPde.getSolution()
295 r=Lsup(length(grad(self.__phi))-1.)
296 vol,vol_old=self.getVolumeOfNegativeDomain(),vol
297 diff=(vol-vol_old)/vol
298 print "iteration :", inf(self.__phi),sup(self.__phi),r,diff
299 # saveVTK("test.%s.xml"%0,l=length(grad(self.__phi,fs)),s=s,phi=self.__phi,v=grad(self.__phi,fs),s2=s2)
300 return
301 #=============================================
302
303 #============
304
305
306 def getVolumeOfNegativeDomain(self):
307 """
308 return the current volume of domain with phi<0.
309 """
310 return integrate((1.-self.__makeInterface(1.))/2.)
311
312
313 def getBoundingBoxOfNegativeDomain(self):
314 """
315 get the height of the region with phi<0
316 """
317 fs=self.__h.getFunctionSpace()
318 mask_phi1=wherePositive(interpolate(self.__phi,fs))
319 mask_phi2=wherePositive(self.__phi)
320 x1=fs.getX()
321 x2=self.__domain.getX()
322 out=[]
323 for i in range(fs.getDim()):
324 x1_i=x1[i]
325 x2_i=x2[i]
326 d=2*(sup(x2_i)-inf(x2_i))
327 offset1=d*mask_phi1
328 offset2=d*mask_phi2
329 out.append((min(inf(x1_i+offset1),inf(x2_i+offset2)),max(sup(x1_i-offset1),sup(x2_i-offset2))))
330 return tuple(out)
331

  ViewVC Help
Powered by ViewVC 1.1.26