/[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 1659 by gross, Fri Jul 18 02:28:13 2008 UTC revision 1661 by gross, Mon Jul 21 22:08:27 2008 UTC
# Line 18  class LevelSet(object): Line 18  class LevelSet(object):
18               xi=x[i]               xi=x[i]
19               diam+=(inf(xi)-sup(xi))**2               diam+=(inf(xi)-sup(xi))**2
20           self.__diam=sqrt(diam)           self.__diam=sqrt(diam)
21           self.__h = Function(self.__domain).getSize()           self.__h = sup(Function(self.__domain).getSize())
22           self.__reinitialization_each=reinitialization_each           self.__reinitialization_each=reinitialization_each
23           self.__reinitialization_steps_max = reinitialization_steps_max           self.__reinitialization_steps_max = reinitialization_steps_max
24           self.__relative_smoothing_width = relative_smoothing_width           self.__relative_smoothing_width = relative_smoothing_width
# Line 28  class LevelSet(object): Line 28  class LevelSet(object):
28    
29           #==================================================           #==================================================
30           self.__fc=TransportPDE(self.__domain,num_equations=1,theta=0.5)           self.__fc=TransportPDE(self.__domain,num_equations=1,theta=0.5)
31             # self.__fc.setReducedOrderOn()
32           self.__fc.setValue(M=Scalar(1.,Function(self.__domain)))           self.__fc.setValue(M=Scalar(1.,Function(self.__domain)))
33           self.__fc.setInitialSolution(phi+self.__diam)           self.__fc.setInitialSolution(phi+self.__diam)
34    
# Line 40  class LevelSet(object): Line 41  class LevelSet(object):
41           # self.__reinitPde.setTolerance(1.e-8)           # self.__reinitPde.setTolerance(1.e-8)
42           # self.__reinitPde.setSolverMethod(preconditioner=LinearPDE.ILU0)           # self.__reinitPde.setSolverMethod(preconditioner=LinearPDE.ILU0)
43           #=======================================           #=======================================
           
   
   
44           self.__updateInterface()           self.__updateInterface()
45             print "phi range:",inf(phi), sup(phi)
46    
47         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    
52       def getDomain(self):       def getDomain(self):
53           return self.__domain           return self.__domain
# Line 53  class LevelSet(object): Line 57  class LevelSet(object):
57           returns a new dt for a given velocity using the courant coundition           returns a new dt for a given velocity using the courant coundition
58           """           """
59           self.velocity=velocity           self.velocity=velocity
60           self.__fc.setValue(C=interpolate(velocity,Function(self.__domain)))           self.__fc.setValue(C=-interpolate(velocity,Function(self.__domain)))
61           dt=self.__fc.getSafeTimeStepSize()           dt=self.__fc.getSafeTimeStepSize()
62           return dt           return dt
63    
# Line 63  class LevelSet(object): Line 67  class LevelSet(object):
67           """           """
68           return self.__phi           return self.__phi
69    
70       def __makeInterface(self,smoothing_width=1.):       def __updateInterface(self):
71             self.__smoothed_char=self.__makeInterface(self.__phi)
72    
73         def __makeInterface(self,phi,smoothing_width=1.):
74           """           """
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           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           and 1 where the level set is 1
77           """           """
          # s=smoothing_width*self.__h*length(grad(self.__phi,self.__h.getFunctionSpace()))  
78           s=smoothing_width*self.__h           s=smoothing_width*self.__h
79           phi_on_h=interpolate(self.__phi,self.__h.getFunctionSpace())                   phi_on_h=interpolate(phi,Function(self.__domain))        
80           mask_neg = whereNonNegative(-s-phi_on_h)           mask_neg = whereNonNegative(-s-phi_on_h)
81           mask_pos = whereNonNegative(phi_on_h-s)           mask_pos = whereNonNegative(phi_on_h-s)
82           mask_interface = 1.-mask_neg-mask_pos           mask_interface = 1.-mask_neg-mask_pos
83           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=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           interface=phi_on_h/s
85           return - mask_neg + mask_pos + mask_interface * interface           return - mask_neg + mask_pos + mask_interface * interface
86    
87         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       def update(self,dt):       def update(self,dt):
99           """           """
100           sets a new velocity and updates the level set fuction           sets a new velocity and updates the level set fuction
101    
102           @param dt: time step forward           @param dt: time step forward
103           """           """
          print inf(self.__phi), sup(self.__phi)  
104           self.__phi=self.__fc.solve(dt, verbose=False)-self.__diam           self.__phi=self.__fc.solve(dt, verbose=False)-self.__diam
          print inf(self.__phi), sup(self.__phi)  
105           self.__update_count += 1           self.__update_count += 1
          # self.__updateInterface()  
106           if self.__update_count%self.__reinitialization_each == 0:           if self.__update_count%self.__reinitialization_each == 0:
107              phi=self.__reinitialise(self.__phi)              self.__phi=self.__reinitialise(self.__phi)
108              self.__fc.setInitialSolution(phi+self.__diam)              self.__fc.setInitialSolution(self.__phi+self.__diam)
109              # self.__updateInterface()           self.__updateInterface()
110    
111       def setTolerance(self,tolerance=1e-3):       def setTolerance(self,tolerance=1e-3):
112           self.__fc.setTolerance(tolerance)           self.__fc.setTolerance(tolerance)
113    
114       def __reinitialise(self,phi):       def __reinitialise(self,phi):
115           print "reintialization started:"           print "reintialization started:"
116           s=self.__makeInterface(1.)           s=self.__makeInterface(phi,1.)
117           g=grad(phi)           g=grad(phi)
118           w = s*g/(length(g)+EPSILON)           w = s*g/(length(g)+EPSILON)
119           dtau = 0.5*inf(self.__h)           dtau = 0.5*inf(Function(self.__domain).getSize())
120           print "step size: dt = ",dtau           print "step size: dt = ",dtau
121           iter =0           iter =0
122             # self.__reinitPde.setValue(q=whereNegative(abs(phi)-self.__h), r=phi)
123             # self.__reinitPde.setValue(r=phi)
124           while (iter<=self.__reinitialization_steps_max):           while (iter<=self.__reinitialization_steps_max):
125                   phi_old=phi                   phi_old=phi
                  print inf(phi+(dtau/2.)*(s-inner(w,grad(phi)))),sup(phi+(dtau/2.)*(s-inner(w,grad(phi))))  
126                   self.__reinitPde.setValue(Y = phi+(dtau/2.)*(s-inner(w,grad(phi))))                   self.__reinitPde.setValue(Y = phi+(dtau/2.)*(s-inner(w,grad(phi))))
127                   phi_half = self.__reinitPde.getSolution()                   phi_half = self.__reinitPde.getSolution()
128                   self.__reinitPde.setValue(Y = phi+dtau*(s-inner(w,grad(phi_half))))                   self.__reinitPde.setValue(Y = phi+dtau*(s-inner(w,grad(phi_half))))
# Line 119  class LevelSet(object): Line 134  class LevelSet(object):
134           print "reintialization completed."           print "reintialization completed."
135           return phi           return phi
136    
137         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    
141    
142       #================ things from here onwards are not used nor tested: ==========================================       #================ things from here onwards are not used nor tested: ==========================================
143            
       
      def __updateInterface(self):  
          self.__smoothed_char=self.__makeInterface(self.__relative_smoothing_width)  
   
144    
145       def getCharacteristicFunction(self):       def getCharacteristicFunction(self):
146           return self.__smoothed_char           return self.__smoothed_char
# Line 201  class LevelSet(object): Line 216  class LevelSet(object):
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)            # 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            self.__phi, previous = self.__reinitPde.getSolution(), self.__phi
219            self.__updateInterface()            # self.__updateInterface()
220    
221            vol,vol_old=self.getVolumeOfNegativeDomain(),vol            vol,vol_old=self.getVolumeOfNegativeDomain(),vol
222            diff=(vol-vol_old)            diff=(vol-vol_old)
# Line 275  class LevelSet(object): Line 290  class LevelSet(object):
290               out.append((min(inf(x1_i+offset1),inf(x2_i+offset2)),max(sup(x1_i-offset1),sup(x2_i-offset2))))               out.append((min(inf(x1_i+offset1),inf(x2_i+offset2)),max(sup(x1_i-offset1),sup(x2_i-offset2))))
291           return tuple(out)           return tuple(out)
292    
      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  

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

  ViewVC Help
Powered by ViewVC 1.1.26