/[escript]/trunk/escript/py_src/levelset.py
ViewVC logotype

Diff of /trunk/escript/py_src/levelset.py

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

revision 1473 by gross, Fri Apr 4 07:49:18 2008 UTC revision 1659 by gross, Fri Jul 18 02:28:13 2008 UTC
# Line 1  Line 1 
1  from esys.escript import *  from esys.escript import *
2  from esys.escript.linearPDEs import LinearPDE  from esys.escript.linearPDEs import LinearPDE, TransportPDE
3  from esys.escript.pdetools import Projector  from esys.escript.pdetools import Projector
4  from esys import finley  from esys import finley
5  import math  import math
# Line 7  import math Line 7  import math
7    
8    
9  class LevelSet(object):  class LevelSet(object):
10       def __init__(self,phi,reinitialization_steps_max=10,reinitialization_each=2,relative_smoothing_width=2., location_fixed_phi=Data()):       def __init__(self,phi,reinitialization_steps_max=10,reinitialization_each=2,relative_smoothing_width=2.):
11           """           """
12           initialize model           initialize model
13           """           """
14           self.domain = phi.getDomain()           self.__domain = phi.getDomain()
15           self.__h = Function(self.domain).getSize()           x=self.__domain.getX()
16           self.__location_fixed_phi=location_fixed_phi           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           self.__relative_smoothing_width = relative_smoothing_width
25           self.__pde = LinearPDE(self.domain)                   self.__phi = phi
26           self.__pde.setReducedOrderOn()           self.__update_count=0
27           self.__pde.setValue(D=1.0,q=self.__location_fixed_phi)           self.velocity = None
28           # self.__pde.setTolerance(1.e-10)  
29           #self.__pde.setSolverMethod(solver=LinearPDE.DIRECT)           #==================================================
30                                 self.__fc=TransportPDE(self.__domain,num_equations=1,theta=0.5)
31           self.__reinitPde = LinearPDE(self.domain,numEquations=1)           self.__fc.setValue(M=Scalar(1.,Function(self.__domain)))
32           # self.__reinitPde.setSolverMethod(solver=LinearPDE.LUMPING)                 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           self.__reinitPde.setReducedOrderOn()           self.__reinitPde.setReducedOrderOn()
          self.__pde.setValue(q=self.__location_fixed_phi)  
38           self.__reinitPde.setSymmetryOn()           self.__reinitPde.setSymmetryOn()
39           self.__reinitPde.setTolerance(1.e-8)           self.__reinitPde.setValue(D=1.)
40             # self.__reinitPde.setTolerance(1.e-8)
41           # self.__reinitPde.setSolverMethod(preconditioner=LinearPDE.ILU0)           # self.__reinitPde.setSolverMethod(preconditioner=LinearPDE.ILU0)
42             #=======================================
43                    
          self.__reinitialization_each=reinitialization_each  
          self.__reinitialization_steps_max = reinitialization_steps_max  
          print self.__reinitialization_steps_max  
          self.__update_count=0  
44    
45    
          self.__phi = phi  
46           self.__updateInterface()           self.__updateInterface()
          self.velocity = None  
47    
48       def getDomain(self):       def getDomain(self):
49           return self.domain           return self.__domain
50        
51             def getTimeStepSize(self,velocity):
52       def __updateInterface(self):           """
53           self.__smoothed_char=self.__makeInterface(self.__relative_smoothing_width)           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       def __makeInterface(self,smoothing_width=1.):       def __makeInterface(self,smoothing_width=1.):
67             """
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           # s=smoothing_width*self.__h*length(grad(self.__phi,self.__h.getFunctionSpace()))           # s=smoothing_width*self.__h*length(grad(self.__phi,self.__h.getFunctionSpace()))
72           s=smoothing_width*self.__h           s=smoothing_width*self.__h
73           phi_on_h=interpolate(self.__phi,self.__h.getFunctionSpace())                   phi_on_h=interpolate(self.__phi,self.__h.getFunctionSpace())        
# Line 58  class LevelSet(object): Line 78  class LevelSet(object):
78           # interface=phi_on_h/s           # interface=phi_on_h/s
79           return - mask_neg + mask_pos + mask_interface * interface           return - mask_neg + mask_pos + mask_interface * interface
80    
      def getCharacteristicFunction(self):  
          return self.__smoothed_char  
   
81       def update(self,dt):       def update(self,dt):
82           """           """
83           sets a new velocity and updates the level set fuction           sets a new velocity and updates the level set fuction
84    
85           @param dt: time step forward           @param dt: time step forward
86           """           """
87           self.__advectPhi(dt)           print inf(self.__phi), sup(self.__phi)
88           self.__updateInterface()           self.__phi=self.__fc.solve(dt, verbose=False)-self.__diam
89           if self.__update_count%self.__reinitialization_each == 0:           print inf(self.__phi), sup(self.__phi)
             self.__reinitialise()  
             self.__updateInterface()  
90           self.__update_count += 1           self.__update_count += 1
91             # self.__updateInterface()
92             if self.__update_count%self.__reinitialization_each == 0:
93                phi=self.__reinitialise(self.__phi)
94                self.__fc.setInitialSolution(phi+self.__diam)
95                # self.__updateInterface()
96    
97         def setTolerance(self,tolerance=1e-3):
98             self.__fc.setTolerance(tolerance)
99    
100         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    
122    
123         #================ 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    
      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")  
134    
      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())))  
135    
136       def RK2(self,L):       def RK2(self,L):
137             k0=L(phi)             k0=L(phi)
# Line 121  class LevelSet(object): Line 152  class LevelSet(object):
152             phi=phi+dt*k1             phi=phi+dt*k1
153                        
154                        
155       def __reinitialise(self):       def __reinitialise_old(self):
156          #=============================================          #=============================================
157          f=0.1          f=0.1
158          f2=0.3          f2=0.3
# Line 163  class LevelSet(object): Line 194  class LevelSet(object):
194            # self.__reinitPde.setValue(Y =self.__phi+dtau/2*s*(1.-len_grad_phi))            # self.__reinitPde.setValue(Y =self.__phi+dtau/2*s*(1.-len_grad_phi))
195            # phi_half=self.__reinitPde.getSolution()            # phi_half=self.__reinitPde.getSolution()
196            # self.__reinitPde.setValue(Y =self.__phi+dtau*s*(1.-length(grad(phi_half,fs))))            # self.__reinitPde.setValue(Y =self.__phi+dtau*s*(1.-length(grad(phi_half,fs))))
197            # 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, Y_reduced=dtau*s*(1.-inner(w,grad(phi_half,ReducedFunction(self.__domain)))))
198            # 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-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)            # 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)            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)
# Line 224  class LevelSet(object): Line 255  class LevelSet(object):
255           """           """
256           return integrate((1.-self.__makeInterface(1.))/2.)           return integrate((1.-self.__makeInterface(1.))/2.)
257    
      def getLevelSetFunction(self):  
          """  
          returns the level set function  
          """  
          return self.__phi  
258                                
259       def getBoundingBoxOfNegativeDomain(self):       def getBoundingBoxOfNegativeDomain(self):
260           """           """
# Line 238  class LevelSet(object): Line 264  class LevelSet(object):
264           mask_phi1=wherePositive(interpolate(self.__phi,fs))           mask_phi1=wherePositive(interpolate(self.__phi,fs))
265           mask_phi2=wherePositive(self.__phi)           mask_phi2=wherePositive(self.__phi)
266           x1=fs.getX()           x1=fs.getX()
267           x2=self.domain.getX()           x2=self.__domain.getX()
268           out=[]           out=[]
269           for i in range(fs.getDim()):           for i in range(fs.getDim()):
270               x1_i=x1[i]               x1_i=x1[i]

Legend:
Removed from v.1473  
changed lines
  Added in v.1659

  ViewVC Help
Powered by ViewVC 1.1.26