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