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

  ViewVC Help
Powered by ViewVC 1.1.26