escript/py_src/levelset.py

from esys.escript import *
from esys.escript.linearPDEs import LinearPDE
from esys.escript.pdetools import Projector
from esys import finley
import math


class LevelSet(object):
     def __init__(self,phi,reinitialization_steps_max=10,reinitialization_each=2,relative_smoothing_width=2., location_fixed_phi=Data()):
         """
         initialize model
         """
         self.domain = phi.getDomain() 
         self.__h = Function(self.domain).getSize()
         self.__location_fixed_phi=location_fixed_phi

         self.__relative_smoothing_width = relative_smoothing_width
         self.__pde = LinearPDE(self.domain)         
         self.__pde.setReducedOrderOn()
         self.__pde.setValue(D=1.0,q=self.__location_fixed_phi)
         # self.__pde.setTolerance(1.e-10)
         #self.__pde.setSolverMethod(solver=LinearPDE.DIRECT)
                     
         self.__reinitPde = LinearPDE(self.domain,numEquations=1)
         # self.__reinitPde.setSolverMethod(solver=LinearPDE.LUMPING)       
         self.__reinitPde.setReducedOrderOn()
         self.__pde.setValue(q=self.__location_fixed_phi)
         self.__reinitPde.setSymmetryOn()
         self.__reinitPde.setTolerance(1.e-8)
         # self.__reinitPde.setSolverMethod(preconditioner=LinearPDE.ILU0)
        
         self.__reinitialization_each=reinitialization_each
         self.__reinitialization_steps_max = reinitialization_steps_max
         print self.__reinitialization_steps_max
         self.__update_count=0


         self.__phi = phi
         self.__updateInterface()
         self.velocity = None

     def getDomain(self):
         return self.domain
     
    
     def __updateInterface(self):
         self.__smoothed_char=self.__makeInterface(self.__relative_smoothing_width)

     def __makeInterface(self,smoothing_width=1.):
         # s=smoothing_width*self.__h*length(grad(self.__phi,self.__h.getFunctionSpace()))
         s=smoothing_width*self.__h
         phi_on_h=interpolate(self.__phi,self.__h.getFunctionSpace())         
         mask_neg = whereNonNegative(-s-phi_on_h)
         mask_pos = whereNonNegative(phi_on_h-s)
         mask_interface = 1.-mask_neg-mask_pos
         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)=
         # interface=phi_on_h/s
         return - mask_neg + mask_pos + mask_interface * interface

     def getCharacteristicFunction(self):
         return self.__smoothed_char

     def update(self,dt):
         """
         sets a new velocity and updates the level set fuction

         @param dt: time step forward
         """
         self.__advectPhi(dt)
         self.__updateInterface()
         if self.__update_count%self.__reinitialization_each == 0:
            self.__reinitialise()
            self.__updateInterface()
         self.__update_count += 1

     def __advectPhi(self,dt):
        """
        advects the level set function 
        """
        fdf=0.01
        grad_phi = grad(self.__phi)   
        len_grad_phi=length(grad_phi)
        f=whereNonNegative(len_grad_phi-fdf)
        s=(1.-f)*len_grad_phi/fdf+f
        
        # self.__pde.setValue(Y=self.__phi-dt/2.*s*inner(self.velocity,grad_phi),r=self.__phi)   
        # phi_half = self.__pde.getSolution()
        # self.__pde.setValue(Y=self.__phi-dt*s*inner(self.velocity,grad(phi_half)),r=self.__phi)   
        # self.__phi = self.__pde.getSolution()

        self.__pde.setValue(Y=s*inner(self.velocity,grad_phi))
        phi_half = self.__phi-dt/2*self.__pde.getSolution()
        self.__pde.setValue(Y=s*inner(self.velocity,grad(phi_half)))
        self.__phi =self.__phi-dt*self.__pde.getSolution()
        print("level set advection done")

     def getTimeStepSize(self,velocity,safety_factor=0.6):
         """
         returns a new dt for a given velocity using the courant coundition
         """
         self.velocity=velocity
         return safety_factor*inf(self.__h/length(interpolate(velocity,self.__h.getFunctionSpace())))

     def RK2(self,L):
           k0=L(phi)
           phi_1=phi+dt*k0
           k1=L(phi_1)
           phi=phi+dt/2*(k0+k1)
     def RK3(self,L):
           k0=L(phi)
           phi_1=phi+dt*L(phi)
           k1=L(phi_1)
           phi_2=phi+dt/4*(k0+k1)
           k2=L(phi_2)
           phi=phi+dt/6*(k0+4*k2+K1)
     def TG(self,L):
           k0=L(phi)
           phi_1=phi+dt/2*k0
           k1=L(phi_1)
           phi=phi+dt*k1
           
           
     def __reinitialise(self):
        #=============================================
        f=0.1
        f2=0.3
        h=self.__h
        fs=h.getFunctionSpace()
        vol=self.getVolumeOfNegativeDomain()
        self.__reinitPde.setValue(D=1.0)
        #============
        grad_phi=grad(self.__phi,fs)
        len_grad_phi=length(grad_phi)
        w=grad_phi/len_grad_phi

        # s=wherePositive(self. __makeInterface(2.))-whereNegative(self. __makeInterface(2.))
        s=self. __makeInterface(2.)
        # s2=whereNegative(len_grad_phi-f2)*(1-(1-len_grad_phi/f2)**2)+(1-whereNegative(len_grad_phi-f2))
        s2=1.
        s=s*s2
        dtau=f*inf(h/abs(s))

        #========================
        diff=1.
        d=1.
        c =0
        #==============================================
        TVD=integrate(length(grad_phi))
        print "initial range ",inf(self.__phi),sup(self.__phi),"error:",Lsup(1.-len_grad_phi),"volume =",vol,TVD
        # saveVTK("test.%s.xml"%c,l=length(grad(self.__phi,fs))-1,s=s,phi=self.__phi)

        dtau=f*inf(h/abs(s))
        while c < self.__reinitialization_steps_max: # and abs(diff) >= 0.01:
          #
          grad_phi=grad(self.__phi,fs)
          len_grad_phi=length(grad_phi)
          #
          # s=self.__makeInterface(1.)
          # s2=whereNegative(len_grad_phi-f2)*(1-(1-len_grad_phi/f2)**2)+(1-whereNegative(len_grad_phi-f2))
          # s=s*s2
          # phi_on_h=interpolate(self.__phi,fs)         
          # self.__reinitPde.setValue(Y =self.__phi+dtau/2*s*(1.-len_grad_phi))
          # phi_half=self.__reinitPde.getSolution()
          # self.__reinitPde.setValue(Y =self.__phi+dtau*s*(1.-length(grad(phi_half,fs))))
          # self.__reinitPde.setValue(Y =self.__phi, Y_reduced=dtau*s*(1.-inner(w,grad(phi_half,ReducedFunction(self.domain)))))
          # self.__reinitPde.setValue(Y =self.__phi+dtau*(s-inner(w,grad_phi)),X=dtau/2*h*(s-inner(w,grad_phi))*w)
          # self.__reinitPde.setValue(Y =self.__phi+dtau*s*(1.-len_grad_phi),X=-dtau/2*h*s**2*grad_phi)
          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)

          # 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)
          self.__phi, previous = self.__reinitPde.getSolution(), self.__phi
          self.__updateInterface()

          vol,vol_old=self.getVolumeOfNegativeDomain(),vol
          diff=(vol-vol_old)
          r=Lsup(length(grad(self.__phi))-1.)
          TVD=integrate(length(grad(self.__phi,fs)))
          print "iteration :", c, "range ",inf(self.__phi),sup(self.__phi),"error :",r,"volume change:",diff,TVD
          # saveVTK("test.%s.xml"%(c+1),l=length(grad(self.__phi,fs)),s=s,phi=self.__phi,v=grad(self.__phi,fs))
          c += 1
        return
        #==============================================
        f2=0.7
        h=self.__h
        fs=h.getFunctionSpace()
        vol=self.getVolumeOfNegativeDomain()
        s=abs(self. __makeInterface(2.))
        grad_phi=grad(self.__phi,fs)
        len_grad_phi=length(grad_phi)
         
        #----
        # aphi_on_h=abs(interpolate(self.__phi,self.__h.getFunctionSpace()))
        # q=2*self.__h
        # mask_neg = whereNegative(aphi_on_h-q)
        # mask_pos = wherePositive(aphi_on_h-q-self.__h)
        # mask_interface = 1.-mask_neg-mask_pos
        # s2=mask_neg + mask_interface * (q+self.__h-aphi_on_h)/self.__h
        #----
        m=whereNonPositive(len_grad_phi-f2)
        s2=m*len_grad_phi/f2+(1.-m)
        #----
        # s2=1.
        #----
        self.__reinitPde.setValue(D=(1-s)/h**2,Y=(1-s)/h**2*self.__phi)
        self.__reinitPde.setValue(A=(s2*len_grad_phi+(1.-s2))*kronecker(3),X=grad_phi)
        # self.__reinitPde.setValue(A=(len_grad_phi-1)/len_grad_phi*kronecker(3))
        # self.__reinitPde.setValue(A=kronecker(3), X=grad_phi/len_grad_phi)
        self.__phi = self.__reinitPde.getSolution()
        r=Lsup(length(grad(self.__phi))-1.)
        vol,vol_old=self.getVolumeOfNegativeDomain(),vol
        diff=(vol-vol_old)/vol
        print "iteration :", inf(self.__phi),sup(self.__phi),r,diff
        # saveVTK("test.%s.xml"%0,l=length(grad(self.__phi,fs)),s=s,phi=self.__phi,v=grad(self.__phi,fs),s2=s2)
        return
        #=============================================

        #============


     def getVolumeOfNegativeDomain(self):
         """
         return the current volume of domain with phi<0.
         """
         return integrate((1.-self.__makeInterface(1.))/2.)

     def getLevelSetFunction(self):
         """
         returns the level set function
         """
         return self.__phi
              
     def getBoundingBoxOfNegativeDomain(self):
         """
         get the height of the region with  phi<0
         """
         fs=self.__h.getFunctionSpace()
         mask_phi1=wherePositive(interpolate(self.__phi,fs))
         mask_phi2=wherePositive(self.__phi)
         x1=fs.getX()
         x2=self.domain.getX()
         out=[]
         for i in range(fs.getDim()):
             x1_i=x1[i]
             x2_i=x2[i]
             d=2*(sup(x2_i)-inf(x2_i))
             offset1=d*mask_phi1
             offset2=d*mask_phi2
             out.append((min(inf(x1_i+offset1),inf(x2_i+offset2)),max(sup(x1_i-offset1),sup(x2_i-offset2))))
         return tuple(out)

     def createParameter(self,value_negative=-1.,value_positive=1):
         out = (value_negative+value_positive)/2. + self.__smoothed_char * ((value_positive-value_negative)/2.)
         return out
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_widthself.__hlength(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.sinner(self.velocity,grad_phi),r=self.__phi)
88	# phi_half = self.__pde.getSolution()
89	# self.__pde.setValue(Y=self.__phi-dtsinner(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+4k2+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/2s(1.-len_grad_phi))
164	# phi_half=self.__reinitPde.getSolution()
165	# self.__reinitPde.setValue(Y =self.__phi+dtaus(1.-length(grad(phi_half,fs))))
166	# self.__reinitPde.setValue(Y =self.__phi, Y_reduced=dtaus(1.-inner(w,grad(phi_half,ReducedFunction(self.domain)))))
167	# self.__reinitPde.setValue(Y =self.__phi+dtau(s-inner(w,grad_phi)),X=dtau/2h(s-inner(w,grad_phi))w)
168	# self.__reinitPde.setValue(Y =self.__phi+dtaus(1.-len_grad_phi),X=-dtau/2hs*2grad_phi)
169	self.__reinitPde.setValue(Y=self.__phi+dtaus(1.-len_grad_phi),X=fdtau/2habs(s)(1.-len_grad_phi)*grad_phi/len_grad_phi)
170
171	# self.__reinitPde.setValue(Y=self.__phi+dtaus(1.-len_grad_phi),X=fdtau/2habs(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)/h2,Y=(1-s)/h2*self.__phi)
206	self.__reinitPde.setValue(A=(s2len_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