/[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 1659 - (show annotations)
Fri Jul 18 02:28:13 2008 UTC (11 years, 1 month ago) by gross
File MIME type: text/x-python
File size: 11074 byte(s)
some first version of a robust 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 = 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=TransportPDE(self.__domain,num_equations=1,theta=0.5)
31 self.__fc.setValue(M=Scalar(1.,Function(self.__domain)))
32 self.__fc.setInitialSolution(phi+self.__diam)
33
34 #=======================================
35 self.__reinitPde = LinearPDE(self.__domain,numEquations=1, numSolutions=1)
36 self.__reinitPde.setSolverMethod(solver=LinearPDE.LUMPING)
37 self.__reinitPde.setReducedOrderOn()
38 self.__reinitPde.setSymmetryOn()
39 self.__reinitPde.setValue(D=1.)
40 # self.__reinitPde.setTolerance(1.e-8)
41 # self.__reinitPde.setSolverMethod(preconditioner=LinearPDE.ILU0)
42 #=======================================
43
44
45
46 self.__updateInterface()
47
48 def getDomain(self):
49 return self.__domain
50
51 def getTimeStepSize(self,velocity):
52 """
53 returns a new dt for a given velocity using the courant coundition
54 """
55 self.velocity=velocity
56 self.__fc.setValue(C=interpolate(velocity,Function(self.__domain)))
57 dt=self.__fc.getSafeTimeStepSize()
58 return dt
59
60 def getLevelSetFunction(self):
61 """
62 returns the level set function
63 """
64 return self.__phi
65
66 def __makeInterface(self,smoothing_width=1.):
67 """
68 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
69 and 1 where the level set is 1
70 """
71 # s=smoothing_width*self.__h*length(grad(self.__phi,self.__h.getFunctionSpace()))
72 s=smoothing_width*self.__h
73 phi_on_h=interpolate(self.__phi,self.__h.getFunctionSpace())
74 mask_neg = whereNonNegative(-s-phi_on_h)
75 mask_pos = whereNonNegative(phi_on_h-s)
76 mask_interface = 1.-mask_neg-mask_pos
77 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)=
78 # interface=phi_on_h/s
79 return - mask_neg + mask_pos + mask_interface * interface
80
81 def update(self,dt):
82 """
83 sets a new velocity and updates the level set fuction
84
85 @param dt: time step forward
86 """
87 print inf(self.__phi), sup(self.__phi)
88 self.__phi=self.__fc.solve(dt, verbose=False)-self.__diam
89 print inf(self.__phi), sup(self.__phi)
90 self.__update_count += 1
91 # self.__updateInterface()
92 if self.__update_count%self.__reinitialization_each == 0:
93 phi=self.__reinitialise(self.__phi)
94 self.__fc.setInitialSolution(phi+self.__diam)
95 # self.__updateInterface()
96
97 def setTolerance(self,tolerance=1e-3):
98 self.__fc.setTolerance(tolerance)
99
100 def __reinitialise(self,phi):
101 print "reintialization started:"
102 s=self.__makeInterface(1.)
103 g=grad(phi)
104 w = s*g/(length(g)+EPSILON)
105 dtau = 0.5*inf(self.__h)
106 print "step size: dt = ",dtau
107 iter =0
108 while (iter<=self.__reinitialization_steps_max):
109 phi_old=phi
110 print inf(phi+(dtau/2.)*(s-inner(w,grad(phi)))),sup(phi+(dtau/2.)*(s-inner(w,grad(phi))))
111 self.__reinitPde.setValue(Y = phi+(dtau/2.)*(s-inner(w,grad(phi))))
112 phi_half = self.__reinitPde.getSolution()
113 self.__reinitPde.setValue(Y = phi+dtau*(s-inner(w,grad(phi_half))))
114 phi = self.__reinitPde.getSolution()
115 change = Lsup(phi-phi_old)/self.__diam
116 print "phi range:",inf(phi), sup(phi)
117 print "iteration :", iter, " error:", change
118 iter +=1
119 print "reintialization completed."
120 return phi
121
122
123 #================ things from here onwards are not used nor tested: ==========================================
124
125
126 def __updateInterface(self):
127 self.__smoothed_char=self.__makeInterface(self.__relative_smoothing_width)
128
129
130 def getCharacteristicFunction(self):
131 return self.__smoothed_char
132
133
134
135
136 def RK2(self,L):
137 k0=L(phi)
138 phi_1=phi+dt*k0
139 k1=L(phi_1)
140 phi=phi+dt/2*(k0+k1)
141 def RK3(self,L):
142 k0=L(phi)
143 phi_1=phi+dt*L(phi)
144 k1=L(phi_1)
145 phi_2=phi+dt/4*(k0+k1)
146 k2=L(phi_2)
147 phi=phi+dt/6*(k0+4*k2+K1)
148 def TG(self,L):
149 k0=L(phi)
150 phi_1=phi+dt/2*k0
151 k1=L(phi_1)
152 phi=phi+dt*k1
153
154
155 def __reinitialise_old(self):
156 #=============================================
157 f=0.1
158 f2=0.3
159 h=self.__h
160 fs=h.getFunctionSpace()
161 vol=self.getVolumeOfNegativeDomain()
162 self.__reinitPde.setValue(D=1.0)
163 #============
164 grad_phi=grad(self.__phi,fs)
165 len_grad_phi=length(grad_phi)
166 w=grad_phi/len_grad_phi
167
168 # s=wherePositive(self. __makeInterface(2.))-whereNegative(self. __makeInterface(2.))
169 s=self. __makeInterface(2.)
170 # s2=whereNegative(len_grad_phi-f2)*(1-(1-len_grad_phi/f2)**2)+(1-whereNegative(len_grad_phi-f2))
171 s2=1.
172 s=s*s2
173 dtau=f*inf(h/abs(s))
174
175 #========================
176 diff=1.
177 d=1.
178 c =0
179 #==============================================
180 TVD=integrate(length(grad_phi))
181 print "initial range ",inf(self.__phi),sup(self.__phi),"error:",Lsup(1.-len_grad_phi),"volume =",vol,TVD
182 # saveVTK("test.%s.xml"%c,l=length(grad(self.__phi,fs))-1,s=s,phi=self.__phi)
183
184 dtau=f*inf(h/abs(s))
185 while c < self.__reinitialization_steps_max: # and abs(diff) >= 0.01:
186 #
187 grad_phi=grad(self.__phi,fs)
188 len_grad_phi=length(grad_phi)
189 #
190 # s=self.__makeInterface(1.)
191 # s2=whereNegative(len_grad_phi-f2)*(1-(1-len_grad_phi/f2)**2)+(1-whereNegative(len_grad_phi-f2))
192 # s=s*s2
193 # phi_on_h=interpolate(self.__phi,fs)
194 # self.__reinitPde.setValue(Y =self.__phi+dtau/2*s*(1.-len_grad_phi))
195 # phi_half=self.__reinitPde.getSolution()
196 # self.__reinitPde.setValue(Y =self.__phi+dtau*s*(1.-length(grad(phi_half,fs))))
197 # self.__reinitPde.setValue(Y =self.__phi, Y_reduced=dtau*s*(1.-inner(w,grad(phi_half,ReducedFunction(self.__domain)))))
198 # self.__reinitPde.setValue(Y =self.__phi+dtau*(s-inner(w,grad_phi)),X=dtau/2*h*(s-inner(w,grad_phi))*w)
199 # self.__reinitPde.setValue(Y =self.__phi+dtau*s*(1.-len_grad_phi),X=-dtau/2*h*s**2*grad_phi)
200 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)
201
202 # 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)
203 self.__phi, previous = self.__reinitPde.getSolution(), self.__phi
204 self.__updateInterface()
205
206 vol,vol_old=self.getVolumeOfNegativeDomain(),vol
207 diff=(vol-vol_old)
208 r=Lsup(length(grad(self.__phi))-1.)
209 TVD=integrate(length(grad(self.__phi,fs)))
210 print "iteration :", c, "range ",inf(self.__phi),sup(self.__phi),"error :",r,"volume change:",diff,TVD
211 # saveVTK("test.%s.xml"%(c+1),l=length(grad(self.__phi,fs)),s=s,phi=self.__phi,v=grad(self.__phi,fs))
212 c += 1
213 return
214 #==============================================
215 f2=0.7
216 h=self.__h
217 fs=h.getFunctionSpace()
218 vol=self.getVolumeOfNegativeDomain()
219 s=abs(self. __makeInterface(2.))
220 grad_phi=grad(self.__phi,fs)
221 len_grad_phi=length(grad_phi)
222
223 #----
224 # aphi_on_h=abs(interpolate(self.__phi,self.__h.getFunctionSpace()))
225 # q=2*self.__h
226 # mask_neg = whereNegative(aphi_on_h-q)
227 # mask_pos = wherePositive(aphi_on_h-q-self.__h)
228 # mask_interface = 1.-mask_neg-mask_pos
229 # s2=mask_neg + mask_interface * (q+self.__h-aphi_on_h)/self.__h
230 #----
231 m=whereNonPositive(len_grad_phi-f2)
232 s2=m*len_grad_phi/f2+(1.-m)
233 #----
234 # s2=1.
235 #----
236 self.__reinitPde.setValue(D=(1-s)/h**2,Y=(1-s)/h**2*self.__phi)
237 self.__reinitPde.setValue(A=(s2*len_grad_phi+(1.-s2))*kronecker(3),X=grad_phi)
238 # self.__reinitPde.setValue(A=(len_grad_phi-1)/len_grad_phi*kronecker(3))
239 # self.__reinitPde.setValue(A=kronecker(3), X=grad_phi/len_grad_phi)
240 self.__phi = self.__reinitPde.getSolution()
241 r=Lsup(length(grad(self.__phi))-1.)
242 vol,vol_old=self.getVolumeOfNegativeDomain(),vol
243 diff=(vol-vol_old)/vol
244 print "iteration :", inf(self.__phi),sup(self.__phi),r,diff
245 # saveVTK("test.%s.xml"%0,l=length(grad(self.__phi,fs)),s=s,phi=self.__phi,v=grad(self.__phi,fs),s2=s2)
246 return
247 #=============================================
248
249 #============
250
251
252 def getVolumeOfNegativeDomain(self):
253 """
254 return the current volume of domain with phi<0.
255 """
256 return integrate((1.-self.__makeInterface(1.))/2.)
257
258
259 def getBoundingBoxOfNegativeDomain(self):
260 """
261 get the height of the region with phi<0
262 """
263 fs=self.__h.getFunctionSpace()
264 mask_phi1=wherePositive(interpolate(self.__phi,fs))
265 mask_phi2=wherePositive(self.__phi)
266 x1=fs.getX()
267 x2=self.__domain.getX()
268 out=[]
269 for i in range(fs.getDim()):
270 x1_i=x1[i]
271 x2_i=x2[i]
272 d=2*(sup(x2_i)-inf(x2_i))
273 offset1=d*mask_phi1
274 offset2=d*mask_phi2
275 out.append((min(inf(x1_i+offset1),inf(x2_i+offset2)),max(sup(x1_i-offset1),sup(x2_i-offset2))))
276 return tuple(out)
277
278 def createParameter(self,value_negative=-1.,value_positive=1):
279 out = (value_negative+value_positive)/2. + self.__smoothed_char * ((value_positive-value_negative)/2.)
280 return out

  ViewVC Help
Powered by ViewVC 1.1.26