/[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 1659 - (hide annotations)
Fri Jul 18 02:28:13 2008 UTC (11 years, 3 months ago) by gross
File MIME type: text/x-python
File size: 11074 byte(s)
some first version of a robust level set
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     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 gross 1473
29 gross 1659 #==================================================
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 gross 1473 self.__reinitPde.setReducedOrderOn()
38     self.__reinitPde.setSymmetryOn()
39 gross 1659 self.__reinitPde.setValue(D=1.)
40     # self.__reinitPde.setTolerance(1.e-8)
41 gross 1473 # self.__reinitPde.setSolverMethod(preconditioner=LinearPDE.ILU0)
42 gross 1659 #=======================================
43 gross 1473
44    
45    
46     self.__updateInterface()
47    
48     def getDomain(self):
49 gross 1659 return self.__domain
50 gross 1473
51 gross 1659 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 gross 1473 def __makeInterface(self,smoothing_width=1.):
67 gross 1659 """
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 gross 1473 # 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 gross 1659 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 gross 1473 if self.__update_count%self.__reinitialization_each == 0:
93 gross 1659 phi=self.__reinitialise(self.__phi)
94     self.__fc.setInitialSolution(phi+self.__diam)
95     # self.__updateInterface()
96 gross 1473
97 gross 1659 def setTolerance(self,tolerance=1e-3):
98     self.__fc.setTolerance(tolerance)
99 gross 1473
100 gross 1659 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 gross 1473
122    
123 gross 1659 #================ 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 gross 1473 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 gross 1659 def __reinitialise_old(self):
156 gross 1473 #=============================================
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 gross 1659 # self.__reinitPde.setValue(Y =self.__phi, Y_reduced=dtau*s*(1.-inner(w,grad(phi_half,ReducedFunction(self.__domain)))))
198 gross 1473 # 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 gross 1659 x2=self.__domain.getX()
268 gross 1473 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