/[escript]/branches/symbolic_from_3470/escript/py_src/levelset.py
ViewVC logotype

Diff of /branches/symbolic_from_3470/escript/py_src/levelset.py

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

revision 3788 by caltinay, Tue Mar 15 04:23:54 2011 UTC revision 3789 by caltinay, Tue Jan 31 04:55:05 2012 UTC
# Line 171  class LevelSet: Line 171  class LevelSet:
171         phi = self.__PDE.getSolution()         phi = self.__PDE.getSolution()
172      else:      else:
173         phi = self.__PDE.getSolution(dt, self.__phi)         phi = self.__PDE.getSolution(dt, self.__phi)
174      print "LevelSet: Advection step done"      print("LevelSet: Advection step done")
175      return phi      return phi
176    
177    #==============================================================================================    #==============================================================================================
# Line 201  class LevelSet: Line 201  class LevelSet:
201        phi = self.__reinitPDE.getSolution()        phi = self.__reinitPDE.getSolution()
202        g=grad(phi)        g=grad(phi)
203        error = Lsup(length(g)*whereNegative(abs(phi.interpolate(g.getFunctionSpace()))-3.0*self.__h))        error = Lsup(length(g)*whereNegative(abs(phi.interpolate(g.getFunctionSpace()))-3.0*self.__h))
204        print "LevelSet:reinitialization: iteration :", iter, " error:", error        print(("LevelSet:reinitialization: iteration :", iter, " error:", error))
205        iter +=1        iter +=1
206      return phi      return phi
207    
# Line 327  class LevelSet2(object): Line 327  class LevelSet2(object):
327              # self.__reinitPde.setSolverMethod(preconditioner=LinearPDE.ILU0)              # self.__reinitPde.setSolverMethod(preconditioner=LinearPDE.ILU0)
328           #=======================================           #=======================================
329           self.__updateInterface()           self.__updateInterface()
330           print "phi range:",inf(phi), sup(phi)           print(("phi range:",inf(phi), sup(phi)))
331    
332       def setFixedRegion(self,mask,contour=0):       def setFixedRegion(self,mask,contour=0):
333           q=whereNonPositive(abs(self.__phi-contour)-self.__h)*mask           q=whereNonPositive(abs(self.__phi-contour)-self.__h)*mask
# Line 364  class LevelSet2(object): Line 364  class LevelSet2(object):
364              self.__fc.setTolerance(tolerance)              self.__fc.setTolerance(tolerance)
365    
366       def __reinitialise(self,phi):       def __reinitialise(self,phi):
367           print "reintialization started:"           print("reintialization started:")
368           s=self.__makeInterface(phi,1.)           s=self.__makeInterface(phi,1.)
369           print "phi range:",inf(phi), sup(phi)           print(("phi range:",inf(phi), sup(phi)))
370           g=grad(phi)           g=grad(phi)
371           w=s*g/(length(g)+EPSILON)           w=s*g/(length(g)+EPSILON)
372           #=======================================================           #=======================================================
# Line 408  class LevelSet2(object): Line 408  class LevelSet2(object):
408               # dtau = 0.5*inf(Function(self.__domain).getSize())               # dtau = 0.5*inf(Function(self.__domain).getSize())
409           else:           else:
410               dtau = 0.5*inf(Function(self.__domain).getSize())               dtau = 0.5*inf(Function(self.__domain).getSize())
411           print "step size: dt = ",dtau,inf(abs(phi.interpolate(Function(self.__domain)))/abs(s-inner(w,grad(phi))))           print(("step size: dt = ",dtau,inf(abs(phi.interpolate(Function(self.__domain)))/abs(s-inner(w,grad(phi))))))
412           iter =0           iter =0
413           # self.__reinitPde.setValue(q=whereNegative(abs(phi)-2*self.__h), r=phi)           # self.__reinitPde.setValue(q=whereNegative(abs(phi)-2*self.__h), r=phi)
414           # self.__reinitPde.setValue(r=phi)           # self.__reinitPde.setValue(r=phi)
# Line 425  class LevelSet2(object): Line 425  class LevelSet2(object):
425                     # self.__reinitPde.setValue(Y = phi+dtau*(s-S),X=dtau**2/2*w*(s-S))                     # self.__reinitPde.setValue(Y = phi+dtau*(s-S),X=dtau**2/2*w*(s-S))
426                     phi = self.__reinitPde.getSolution()                     phi = self.__reinitPde.getSolution()
427                   change = Lsup(phi-phi_old)/self.__diam                   change = Lsup(phi-phi_old)/self.__diam
428                   print "phi range:",inf(phi), sup(phi)                   print(("phi range:",inf(phi), sup(phi)))
429                   print "iteration :", iter, " change:", change                   print(("iteration :", iter, " change:", change))
430                   iter +=1                   iter +=1
431           print "reintialization completed."           print("reintialization completed.")
432           return phi           return phi
433    
434       def createParameter(self,value_negative=-1.,value_positive=1):       def createParameter(self,value_negative=-1.,value_positive=1):
# Line 490  class LevelSet2(object): Line 490  class LevelSet2(object):
490          c =0          c =0
491          #==============================================          #==============================================
492          TVD=integrate(length(grad_phi))          TVD=integrate(length(grad_phi))
493          print "initial range ",inf(self.__phi),sup(self.__phi),"error:",Lsup(1.-len_grad_phi),"volume =",vol,TVD          print(("initial range ",inf(self.__phi),sup(self.__phi),"error:",Lsup(1.-len_grad_phi),"volume =",vol,TVD))
494          # saveVTK("test.%s.vtu"%c,l=length(grad(self.__phi,fs))-1,s=s,phi=self.__phi)          # saveVTK("test.%s.vtu"%c,l=length(grad(self.__phi,fs))-1,s=s,phi=self.__phi)
495    
496          dtau=f*inf(h/abs(s))          dtau=f*inf(h/abs(s))
# Line 519  class LevelSet2(object): Line 519  class LevelSet2(object):
519            diff=(vol-vol_old)            diff=(vol-vol_old)
520            r=Lsup(length(grad(self.__phi))-1.)            r=Lsup(length(grad(self.__phi))-1.)
521            TVD=integrate(length(grad(self.__phi,fs)))            TVD=integrate(length(grad(self.__phi,fs)))
522            print "iteration :", c, "range ",inf(self.__phi),sup(self.__phi),"error :",r,"volume change:",diff,TVD            print(("iteration :", c, "range ",inf(self.__phi),sup(self.__phi),"error :",r,"volume change:",diff,TVD))
523            # saveVTK("test.%s.vtu"%(c+1),l=length(grad(self.__phi,fs)),s=s,phi=self.__phi,v=grad(self.__phi,fs))            # saveVTK("test.%s.vtu"%(c+1),l=length(grad(self.__phi,fs)),s=s,phi=self.__phi,v=grad(self.__phi,fs))
524            c += 1            c += 1
525          return          return
# Line 553  class LevelSet2(object): Line 553  class LevelSet2(object):
553          r=Lsup(length(grad(self.__phi))-1.)          r=Lsup(length(grad(self.__phi))-1.)
554          vol,vol_old=self.getVolumeOfNegativeDomain(),vol          vol,vol_old=self.getVolumeOfNegativeDomain(),vol
555          diff=(vol-vol_old)/vol          diff=(vol-vol_old)/vol
556          print "iteration :", inf(self.__phi),sup(self.__phi),r,diff          print(("iteration :", inf(self.__phi),sup(self.__phi),r,diff))
557          # saveVTK("test.%s.vtu"%0,l=length(grad(self.__phi,fs)),s=s,phi=self.__phi,v=grad(self.__phi,fs),s2=s2)          # saveVTK("test.%s.vtu"%0,l=length(grad(self.__phi,fs)),s=s,phi=self.__phi,v=grad(self.__phi,fs),s2=s2)
558          return          return
559          #=============================================          #=============================================

Legend:
Removed from v.3788  
changed lines
  Added in v.3789

  ViewVC Help
Powered by ViewVC 1.1.26