/[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 1473 - (show annotations)
Fri Apr 4 07:49:18 2008 UTC (11 years, 6 months ago) by gross
File MIME type: text/x-python
File size: 10134 byte(s)
initial checkin for level set
1 from esys.escript import *
2 from esys.escript.linearPDEs import LinearPDE
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., location_fixed_phi=Data()):
11 """
12 initialize model
13 """
14 self.domain = phi.getDomain()
15 self.__h = Function(self.domain).getSize()
16 self.__location_fixed_phi=location_fixed_phi
17
18 self.__relative_smoothing_width = relative_smoothing_width
19 self.__pde = LinearPDE(self.domain)
20 self.__pde.setReducedOrderOn()
21 self.__pde.setValue(D=1.0,q=self.__location_fixed_phi)
22 # self.__pde.setTolerance(1.e-10)
23 #self.__pde.setSolverMethod(solver=LinearPDE.DIRECT)
24
25 self.__reinitPde = LinearPDE(self.domain,numEquations=1)
26 # self.__reinitPde.setSolverMethod(solver=LinearPDE.LUMPING)
27 self.__reinitPde.setReducedOrderOn()
28 self.__pde.setValue(q=self.__location_fixed_phi)
29 self.__reinitPde.setSymmetryOn()
30 self.__reinitPde.setTolerance(1.e-8)
31 # self.__reinitPde.setSolverMethod(preconditioner=LinearPDE.ILU0)
32
33 self.__reinitialization_each=reinitialization_each
34 self.__reinitialization_steps_max = reinitialization_steps_max
35 print self.__reinitialization_steps_max
36 self.__update_count=0
37
38
39 self.__phi = phi
40 self.__updateInterface()
41 self.velocity = None
42
43 def getDomain(self):
44 return self.domain
45
46
47 def __updateInterface(self):
48 self.__smoothed_char=self.__makeInterface(self.__relative_smoothing_width)
49
50 def __makeInterface(self,smoothing_width=1.):
51 # s=smoothing_width*self.__h*length(grad(self.__phi,self.__h.getFunctionSpace()))
52 s=smoothing_width*self.__h
53 phi_on_h=interpolate(self.__phi,self.__h.getFunctionSpace())
54 mask_neg = whereNonNegative(-s-phi_on_h)
55 mask_pos = whereNonNegative(phi_on_h-s)
56 mask_interface = 1.-mask_neg-mask_pos
57 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)=
58 # interface=phi_on_h/s
59 return - mask_neg + mask_pos + mask_interface * interface
60
61 def getCharacteristicFunction(self):
62 return self.__smoothed_char
63
64 def update(self,dt):
65 """
66 sets a new velocity and updates the level set fuction
67
68 @param dt: time step forward
69 """
70 self.__advectPhi(dt)
71 self.__updateInterface()
72 if self.__update_count%self.__reinitialization_each == 0:
73 self.__reinitialise()
74 self.__updateInterface()
75 self.__update_count += 1
76
77 def __advectPhi(self,dt):
78 """
79 advects the level set function
80 """
81 fdf=0.01
82 grad_phi = grad(self.__phi)
83 len_grad_phi=length(grad_phi)
84 f=whereNonNegative(len_grad_phi-fdf)
85 s=(1.-f)*len_grad_phi/fdf+f
86
87 # self.__pde.setValue(Y=self.__phi-dt/2.*s*inner(self.velocity,grad_phi),r=self.__phi)
88 # phi_half = self.__pde.getSolution()
89 # self.__pde.setValue(Y=self.__phi-dt*s*inner(self.velocity,grad(phi_half)),r=self.__phi)
90 # self.__phi = self.__pde.getSolution()
91
92 self.__pde.setValue(Y=s*inner(self.velocity,grad_phi))
93 phi_half = self.__phi-dt/2*self.__pde.getSolution()
94 self.__pde.setValue(Y=s*inner(self.velocity,grad(phi_half)))
95 self.__phi =self.__phi-dt*self.__pde.getSolution()
96 print("level set advection done")
97
98 def getTimeStepSize(self,velocity,safety_factor=0.6):
99 """
100 returns a new dt for a given velocity using the courant coundition
101 """
102 self.velocity=velocity
103 return safety_factor*inf(self.__h/length(interpolate(velocity,self.__h.getFunctionSpace())))
104
105 def RK2(self,L):
106 k0=L(phi)
107 phi_1=phi+dt*k0
108 k1=L(phi_1)
109 phi=phi+dt/2*(k0+k1)
110 def RK3(self,L):
111 k0=L(phi)
112 phi_1=phi+dt*L(phi)
113 k1=L(phi_1)
114 phi_2=phi+dt/4*(k0+k1)
115 k2=L(phi_2)
116 phi=phi+dt/6*(k0+4*k2+K1)
117 def TG(self,L):
118 k0=L(phi)
119 phi_1=phi+dt/2*k0
120 k1=L(phi_1)
121 phi=phi+dt*k1
122
123
124 def __reinitialise(self):
125 #=============================================
126 f=0.1
127 f2=0.3
128 h=self.__h
129 fs=h.getFunctionSpace()
130 vol=self.getVolumeOfNegativeDomain()
131 self.__reinitPde.setValue(D=1.0)
132 #============
133 grad_phi=grad(self.__phi,fs)
134 len_grad_phi=length(grad_phi)
135 w=grad_phi/len_grad_phi
136
137 # s=wherePositive(self. __makeInterface(2.))-whereNegative(self. __makeInterface(2.))
138 s=self. __makeInterface(2.)
139 # s2=whereNegative(len_grad_phi-f2)*(1-(1-len_grad_phi/f2)**2)+(1-whereNegative(len_grad_phi-f2))
140 s2=1.
141 s=s*s2
142 dtau=f*inf(h/abs(s))
143
144 #========================
145 diff=1.
146 d=1.
147 c =0
148 #==============================================
149 TVD=integrate(length(grad_phi))
150 print "initial range ",inf(self.__phi),sup(self.__phi),"error:",Lsup(1.-len_grad_phi),"volume =",vol,TVD
151 # saveVTK("test.%s.xml"%c,l=length(grad(self.__phi,fs))-1,s=s,phi=self.__phi)
152
153 dtau=f*inf(h/abs(s))
154 while c < self.__reinitialization_steps_max: # and abs(diff) >= 0.01:
155 #
156 grad_phi=grad(self.__phi,fs)
157 len_grad_phi=length(grad_phi)
158 #
159 # s=self.__makeInterface(1.)
160 # s2=whereNegative(len_grad_phi-f2)*(1-(1-len_grad_phi/f2)**2)+(1-whereNegative(len_grad_phi-f2))
161 # s=s*s2
162 # phi_on_h=interpolate(self.__phi,fs)
163 # self.__reinitPde.setValue(Y =self.__phi+dtau/2*s*(1.-len_grad_phi))
164 # phi_half=self.__reinitPde.getSolution()
165 # self.__reinitPde.setValue(Y =self.__phi+dtau*s*(1.-length(grad(phi_half,fs))))
166 # self.__reinitPde.setValue(Y =self.__phi, Y_reduced=dtau*s*(1.-inner(w,grad(phi_half,ReducedFunction(self.domain)))))
167 # self.__reinitPde.setValue(Y =self.__phi+dtau*(s-inner(w,grad_phi)),X=dtau/2*h*(s-inner(w,grad_phi))*w)
168 # self.__reinitPde.setValue(Y =self.__phi+dtau*s*(1.-len_grad_phi),X=-dtau/2*h*s**2*grad_phi)
169 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)
170
171 # 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)
172 self.__phi, previous = self.__reinitPde.getSolution(), self.__phi
173 self.__updateInterface()
174
175 vol,vol_old=self.getVolumeOfNegativeDomain(),vol
176 diff=(vol-vol_old)
177 r=Lsup(length(grad(self.__phi))-1.)
178 TVD=integrate(length(grad(self.__phi,fs)))
179 print "iteration :", c, "range ",inf(self.__phi),sup(self.__phi),"error :",r,"volume change:",diff,TVD
180 # saveVTK("test.%s.xml"%(c+1),l=length(grad(self.__phi,fs)),s=s,phi=self.__phi,v=grad(self.__phi,fs))
181 c += 1
182 return
183 #==============================================
184 f2=0.7
185 h=self.__h
186 fs=h.getFunctionSpace()
187 vol=self.getVolumeOfNegativeDomain()
188 s=abs(self. __makeInterface(2.))
189 grad_phi=grad(self.__phi,fs)
190 len_grad_phi=length(grad_phi)
191
192 #----
193 # aphi_on_h=abs(interpolate(self.__phi,self.__h.getFunctionSpace()))
194 # q=2*self.__h
195 # mask_neg = whereNegative(aphi_on_h-q)
196 # mask_pos = wherePositive(aphi_on_h-q-self.__h)
197 # mask_interface = 1.-mask_neg-mask_pos
198 # s2=mask_neg + mask_interface * (q+self.__h-aphi_on_h)/self.__h
199 #----
200 m=whereNonPositive(len_grad_phi-f2)
201 s2=m*len_grad_phi/f2+(1.-m)
202 #----
203 # s2=1.
204 #----
205 self.__reinitPde.setValue(D=(1-s)/h**2,Y=(1-s)/h**2*self.__phi)
206 self.__reinitPde.setValue(A=(s2*len_grad_phi+(1.-s2))*kronecker(3),X=grad_phi)
207 # self.__reinitPde.setValue(A=(len_grad_phi-1)/len_grad_phi*kronecker(3))
208 # self.__reinitPde.setValue(A=kronecker(3), X=grad_phi/len_grad_phi)
209 self.__phi = self.__reinitPde.getSolution()
210 r=Lsup(length(grad(self.__phi))-1.)
211 vol,vol_old=self.getVolumeOfNegativeDomain(),vol
212 diff=(vol-vol_old)/vol
213 print "iteration :", inf(self.__phi),sup(self.__phi),r,diff
214 # saveVTK("test.%s.xml"%0,l=length(grad(self.__phi,fs)),s=s,phi=self.__phi,v=grad(self.__phi,fs),s2=s2)
215 return
216 #=============================================
217
218 #============
219
220
221 def getVolumeOfNegativeDomain(self):
222 """
223 return the current volume of domain with phi<0.
224 """
225 return integrate((1.-self.__makeInterface(1.))/2.)
226
227 def getLevelSetFunction(self):
228 """
229 returns the level set function
230 """
231 return self.__phi
232
233 def getBoundingBoxOfNegativeDomain(self):
234 """
235 get the height of the region with phi<0
236 """
237 fs=self.__h.getFunctionSpace()
238 mask_phi1=wherePositive(interpolate(self.__phi,fs))
239 mask_phi2=wherePositive(self.__phi)
240 x1=fs.getX()
241 x2=self.domain.getX()
242 out=[]
243 for i in range(fs.getDim()):
244 x1_i=x1[i]
245 x2_i=x2[i]
246 d=2*(sup(x2_i)-inf(x2_i))
247 offset1=d*mask_phi1
248 offset2=d*mask_phi2
249 out.append((min(inf(x1_i+offset1),inf(x2_i+offset2)),max(sup(x1_i-offset1),sup(x2_i-offset2))))
250 return tuple(out)
251
252 def createParameter(self,value_negative=-1.,value_positive=1):
253 out = (value_negative+value_positive)/2. + self.__smoothed_char * ((value_positive-value_negative)/2.)
254 return out

  ViewVC Help
Powered by ViewVC 1.1.26