27 |
self.velocity = None |
self.velocity = None |
28 |
|
|
29 |
#================================================== |
#================================================== |
30 |
self.__fc=TransportPDE(self.__domain,num_equations=1,theta=0.5) |
self.__FC=True |
31 |
# self.__fc.setReducedOrderOn() |
if self.__FC: |
32 |
self.__fc.setValue(M=Scalar(1.,Function(self.__domain))) |
self.__fc=TransportPDE(self.__domain,num_equations=1,theta=0.5) |
33 |
self.__fc.setInitialSolution(phi+self.__diam) |
# self.__fc.setReducedOrderOn() |
34 |
|
self.__fc.setValue(M=Scalar(1.,Function(self.__domain))) |
35 |
|
self.__fc.setInitialSolution(phi+self.__diam) |
36 |
|
else: |
37 |
|
self.__fcpde = LinearPDE(self.__domain,numEquations=1, numSolutions=1) |
38 |
|
self.__fcpde.setSolverMethod(solver=LinearPDE.LUMPING) |
39 |
|
self.__fcpde.setReducedOrderOn() |
40 |
|
self.__fcpde.setSymmetryOn() |
41 |
|
self.__fcpde.setValue(D=1.) |
42 |
|
|
43 |
#======================================= |
#======================================= |
44 |
self.__reinitPde = LinearPDE(self.__domain,numEquations=1, numSolutions=1) |
self.__reinitFC=False |
45 |
self.__reinitPde.setSolverMethod(solver=LinearPDE.LUMPING) |
if self.__reinitFC: |
46 |
self.__reinitPde.setReducedOrderOn() |
self.__reinitfc=TransportPDE(self.__domain,num_equations=1,theta=1.) |
47 |
self.__reinitPde.setSymmetryOn() |
self.__reinitfc.setValue(M=Scalar(1.,Function(self.__domain))) |
48 |
self.__reinitPde.setValue(D=1.) |
self.__reinitfc.setTolerance(1.e-5) |
49 |
# self.__reinitPde.setTolerance(1.e-8) |
else: |
50 |
# self.__reinitPde.setSolverMethod(preconditioner=LinearPDE.ILU0) |
self.__reinitPde = LinearPDE(self.__domain,numEquations=1, numSolutions=1) |
51 |
|
# self.__reinitPde.setSolverMethod(solver=LinearPDE.LUMPING) |
52 |
|
self.__reinitPde.setReducedOrderOn() |
53 |
|
self.__reinitPde.setSymmetryOn() |
54 |
|
self.__reinitPde.setValue(D=1.) |
55 |
|
# self.__reinitPde.setTolerance(1.e-8) |
56 |
|
# self.__reinitPde.setSolverMethod(preconditioner=LinearPDE.ILU0) |
57 |
#======================================= |
#======================================= |
58 |
self.__updateInterface() |
self.__updateInterface() |
59 |
print "phi range:",inf(phi), sup(phi) |
print "phi range:",inf(phi), sup(phi) |
63 |
# self.__fc.setValue(q=q) |
# self.__fc.setValue(q=q) |
64 |
self.__reinitPde.setValue(q=q) |
self.__reinitPde.setValue(q=q) |
65 |
|
|
66 |
|
def getH(self): |
67 |
|
return self.__h |
68 |
def getDomain(self): |
def getDomain(self): |
69 |
return self.__domain |
return self.__domain |
70 |
|
|
73 |
returns a new dt for a given velocity using the courant coundition |
returns a new dt for a given velocity using the courant coundition |
74 |
""" |
""" |
75 |
self.velocity=velocity |
self.velocity=velocity |
76 |
self.__fc.setValue(C=-interpolate(velocity,Function(self.__domain))) |
if self.__FC: |
77 |
dt=self.__fc.getSafeTimeStepSize() |
self.__fc.setValue(C=-interpolate(velocity,Function(self.__domain))) |
78 |
|
dt=self.__fc.getSafeTimeStepSize() |
79 |
|
else: |
80 |
|
dt=0.5*self.__h/sup(length(velocity)) |
81 |
return dt |
return dt |
82 |
|
|
83 |
def getLevelSetFunction(self): |
def getLevelSetFunction(self): |
120 |
|
|
121 |
@param dt: time step forward |
@param dt: time step forward |
122 |
""" |
""" |
123 |
self.__phi=self.__fc.solve(dt, verbose=False)-self.__diam |
if self.__FC: |
124 |
|
self.__phi=self.__fc.solve(dt, verbose=False)-self.__diam |
125 |
|
else: |
126 |
|
self.__fcpde.setValue(Y = self.__phi-(dt/2.)*inner(self.velocity,grad(self.__phi))) |
127 |
|
phi_half = self.__fcpde.getSolution() |
128 |
|
self.__fcpde.setValue(Y = self.__phi-dt*inner(self.velocity,grad(phi_half))) |
129 |
|
self.__phi= self.__fcpde.getSolution() |
130 |
self.__update_count += 1 |
self.__update_count += 1 |
131 |
if self.__update_count%self.__reinitialization_each == 0: |
if self.__update_count%self.__reinitialization_each == 0: |
132 |
self.__phi=self.__reinitialise(self.__phi) |
self.__phi=self.__reinitialise(self.__phi) |
133 |
self.__fc.setInitialSolution(self.__phi+self.__diam) |
if self.__FC: |
134 |
|
self.__fc.setInitialSolution(self.__phi+self.__diam) |
135 |
self.__updateInterface() |
self.__updateInterface() |
136 |
|
|
137 |
def setTolerance(self,tolerance=1e-3): |
def setTolerance(self,tolerance=1e-3): |
138 |
self.__fc.setTolerance(tolerance) |
if self.__FC: |
139 |
|
self.__fc.setTolerance(tolerance) |
140 |
|
|
141 |
def __reinitialise(self,phi): |
def __reinitialise(self,phi): |
142 |
print "reintialization started:" |
print "reintialization started:" |
143 |
s=self.__makeInterface(phi,1.) |
s=self.__makeInterface(phi,1.5) |
144 |
g=grad(phi) |
g=grad(phi) |
145 |
w = s*g/(length(g)+EPSILON) |
w = s*g/(length(g)+EPSILON) |
146 |
dtau = 0.5*inf(Function(self.__domain).getSize()) |
if self.__reinitFC: |
147 |
print "step size: dt = ",dtau |
self.__reinitfc.setValue(C=-w,Y=s) |
148 |
|
self.__reinitfc.setInitialSolution(phi+self.__diam) |
149 |
|
# dtau=self.__reinitfc.getSafeTimeStepSize() |
150 |
|
dtau = 0.5*inf(Function(self.__domain).getSize()) |
151 |
|
else: |
152 |
|
dtau = 0.5*inf(Function(self.__domain).getSize()) |
153 |
|
print "step size: dt = ",dtau,inf(abs(phi.interpolate(Function(self.__domain)))/abs(s-inner(w,grad(phi)))) |
154 |
iter =0 |
iter =0 |
155 |
# self.__reinitPde.setValue(q=whereNegative(abs(phi)-self.__h), r=phi) |
# self.__reinitPde.setValue(q=whereNegative(abs(phi)-2*self.__h), r=phi) |
156 |
# self.__reinitPde.setValue(r=phi) |
# self.__reinitPde.setValue(r=phi) |
157 |
while (iter<=self.__reinitialization_steps_max): |
while (iter<=self.__reinitialization_steps_max): |
158 |
phi_old=phi |
phi_old=phi |
159 |
self.__reinitPde.setValue(Y = phi+(dtau/2.)*(s-inner(w,grad(phi)))) |
if self.__reinitFC: |
160 |
phi_half = self.__reinitPde.getSolution() |
phi = self.__reinitfc.solve(dtau)-self.__diam |
161 |
self.__reinitPde.setValue(Y = phi+dtau*(s-inner(w,grad(phi_half)))) |
else: |
162 |
phi = self.__reinitPde.getSolution() |
self.__reinitPde.setValue(Y = phi+(dtau/2.)*(s-inner(w,grad(phi)))) |
163 |
|
phi_half = self.__reinitPde.getSolution() |
164 |
|
self.__reinitPde.setValue(Y = phi+dtau*(s-inner(w,grad(phi_half)))) |
165 |
|
# g=grad(phi) |
166 |
|
# S=inner(w,grad(phi)) |
167 |
|
# self.__reinitPde.setValue(Y = phi+dtau*(s-S),X=dtau**2/2*w*(s-S)) |
168 |
|
phi = self.__reinitPde.getSolution() |
169 |
change = Lsup(phi-phi_old)/self.__diam |
change = Lsup(phi-phi_old)/self.__diam |
170 |
print "phi range:",inf(phi), sup(phi) |
print "phi range:",inf(phi), sup(phi) |
171 |
print "iteration :", iter, " error:", change |
print "iteration :", iter, " error:", change |