/[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 1473 - (hide 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 gross 1473 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