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

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

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

revision 1780 by jfenwick, Mon Aug 11 03:33:40 2008 UTC revision 1781 by jfenwick, Thu Sep 11 05:03:14 2008 UTC
# Line 43  class LevelSet(object): Line 43  class LevelSet(object):
43           #=======================================           #=======================================
44           self.__reinitFC=False           self.__reinitFC=False
45           if self.__reinitFC:           if self.__reinitFC:
46              self.__reinitfc=TransportPDE(self.__domain,num_equations=1,theta=1.)              self.__reinitfc=TransportPDE(self.__domain,num_equations=1,theta=1.0)
47              self.__reinitfc.setValue(M=Scalar(1.,Function(self.__domain)))              self.__reinitfc.setValue(M=Scalar(1.,Function(self.__domain)))
48              self.__reinitfc.setTolerance(1.e-5)              self.__reinitfc.setTolerance(1.e-5)
49           else:           else:
# Line 52  class LevelSet(object): Line 52  class LevelSet(object):
52              self.__reinitPde.setReducedOrderOn()              self.__reinitPde.setReducedOrderOn()
53              self.__reinitPde.setSymmetryOn()              self.__reinitPde.setSymmetryOn()
54              self.__reinitPde.setValue(D=1.)              self.__reinitPde.setValue(D=1.)
55              # self.__reinitPde.setTolerance(1.e-8)              self.__reinitPde.setTolerance(1.e-2)
56              # self.__reinitPde.setSolverMethod(preconditioner=LinearPDE.ILU0)              # self.__reinitPde.setSolverMethod(preconditioner=LinearPDE.ILU0)
57           #=======================================           #=======================================
58           self.__updateInterface()           self.__updateInterface()
# Line 140  class LevelSet(object): Line 140  class LevelSet(object):
140    
141       def __reinitialise(self,phi):       def __reinitialise(self,phi):
142           print "reintialization started:"           print "reintialization started:"
143           s=self.__makeInterface(phi,1.5)           s=self.__makeInterface(phi,1.)
144             print "phi range:",inf(phi), sup(phi)
145           g=grad(phi)           g=grad(phi)
146           w = s*g/(length(g)+EPSILON)           w=s*g/(length(g)+EPSILON)
147             #=======================================================
148             # # positive part:
149             # phi_p=wherePositive(phi)*phi
150             # self.__reinitfc.setInitialSolution(phi_p)
151             # self.__reinitfc.setValue(C=-wherePositive(s)*w,Y=wherePositive(s)*s,q=whereNonPositive(phi))
152             # dtau=self.__h
153             # print "step size: dt (pos)= ",dtau
154             # print "phi_p range:",inf(phi_p), sup(phi_p)
155             # iter=0
156             # while (iter<=self.__reinitialization_steps_max):
157             # phi_p=self.__reinitfc.solve(dtau)
158             # print "phi_p range:",inf(phi_p), sup(phi_p)
159             # iter+=1
160             # # negative part:
161             # phi_n=-whereNegative(phi)*phi
162             # self.__reinitfc.setInitialSolution(phi_n)
163             # self.__reinitfc.setValue(C=-whereNegative(s)*w,Y=-whereNegative(s)*s,q=whereNonNegative(phi))
164             # # dtau=self.__reinitfc.getSafeTimeStepSize()
165             # dtau=self.__h
166             # print "step size: dt (neg)= ",dtau
167             # print "phi_n range:",inf(phi_n), sup(phi_n)
168             # iter=0
169             # while (iter<=self.__reinitialization_steps_max):
170             # phi_n=self.__reinitfc.solve(dtau)
171             # print "phi_n range:",inf(phi_n), sup(phi_n)
172             # iter+=1
173             # phi=phi_p-phi_n
174             # print "phi range:",inf(phi), sup(phi)
175             # print "reintialization completed."
176             # return phi
177    
178             #=======================================================
179           if self.__reinitFC:           if self.__reinitFC:
180               self.__reinitfc.setValue(C=-w,Y=s)               self.__reinitfc.setValue(C=-w,Y=s)
181               self.__reinitfc.setInitialSolution(phi+self.__diam)               self.__reinitfc.setInitialSolution(phi+self.__diam)
182               # dtau=self.__reinitfc.getSafeTimeStepSize()               dtau=self.__reinitfc.getSafeTimeStepSize()
183               dtau = 0.5*inf(Function(self.__domain).getSize())               # dtau = 0.5*inf(Function(self.__domain).getSize())
184           else:           else:
185               dtau = 0.5*inf(Function(self.__domain).getSize())               dtau = 0.5*inf(Function(self.__domain).getSize())
186           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))))
# Line 168  class LevelSet(object): Line 201  class LevelSet(object):
201                     phi = self.__reinitPde.getSolution()                     phi = self.__reinitPde.getSolution()
202                   change = Lsup(phi-phi_old)/self.__diam                   change = Lsup(phi-phi_old)/self.__diam
203                   print "phi range:",inf(phi), sup(phi)                   print "phi range:",inf(phi), sup(phi)
204                   print "iteration :", iter, " error:", change                   print "iteration :", iter, " change:", change
205                   iter +=1                   iter +=1
206           print "reintialization completed."           print "reintialization completed."
207           return phi           return phi

Legend:
Removed from v.1780  
changed lines
  Added in v.1781

  ViewVC Help
Powered by ViewVC 1.1.26