| 1 |
|
| 2 |
######################################################## |
| 3 |
# |
| 4 |
# Copyright (c) 2003-2008 by University of Queensland |
| 5 |
# Earth Systems Science Computational Center (ESSCC) |
| 6 |
# http://www.uq.edu.au/esscc |
| 7 |
# |
| 8 |
# Primary Business: Queensland, Australia |
| 9 |
# Licensed under the Open Software License version 3.0 |
| 10 |
# http://www.opensource.org/licenses/osl-3.0.php |
| 11 |
# |
| 12 |
######################################################## |
| 13 |
|
| 14 |
__copyright__="""Copyright (c) 2003-2008 by University of Queensland |
| 15 |
Earth Systems Science Computational Center (ESSCC) |
| 16 |
http://www.uq.edu.au/esscc |
| 17 |
Primary Business: Queensland, Australia""" |
| 18 |
__license__="""Licensed under the Open Software License version 3.0 |
| 19 |
http://www.opensource.org/licenses/osl-3.0.php""" |
| 20 |
__url__="http://www.uq.edu.au/esscc/escript-finley" |
| 21 |
|
| 22 |
from esys.escript import * |
| 23 |
from esys.escript.linearPDEs import LinearPDE, TransportPDE |
| 24 |
from esys.escript.pdetools import Projector |
| 25 |
from esys import finley |
| 26 |
import math |
| 27 |
|
| 28 |
|
| 29 |
|
| 30 |
from esys.escript import * |
| 31 |
from esys.finley import finley |
| 32 |
from esys.escript.linearPDEs import LinearPDE |
| 33 |
from esys.escript.pdetools import Projector |
| 34 |
import sys |
| 35 |
|
| 36 |
|
| 37 |
class LevelSet: |
| 38 |
""" |
| 39 |
The level set method tracking an interface defined by the zero contour of the |
| 40 |
level set function phi. |
| 41 |
|
| 42 |
It is assumed that phi(x)<0 defines the volume of interest. |
| 43 |
|
| 44 |
""" |
| 45 |
def __init__(self,domain,phi,reinit_max=10,reinit_each=2,smooth=2.): |
| 46 |
""" |
| 47 |
Sets up the level set method. |
| 48 |
|
| 49 |
@param domain: the domain where the level set is used |
| 50 |
@param phi: the initial level set function |
| 51 |
@param reinit_max: maximum number of reinitialization steps |
| 52 |
@param reinit_each: C{phi} is reinitialized every C{reinit_each} step |
| 53 |
@param smooth: smoothing width |
| 54 |
""" |
| 55 |
self.__domain = domain |
| 56 |
self.__phi = phi |
| 57 |
self.__reinit_max = reinit_max |
| 58 |
self.__reinit_each = reinit_each |
| 59 |
self.__PDE = LinearPDE(domain) |
| 60 |
self.__PDE.setReducedOrderOn() |
| 61 |
self.__PDE.setValue(D=1.0) |
| 62 |
self.__PDE.setSolverMethod(solver=LinearPDE.PCG) |
| 63 |
self.__reinitPDE = LinearPDE(domain, numEquations=1) |
| 64 |
self.__reinitPDE.setReducedOrderOn() |
| 65 |
self.__reinitPDE.setValue(D=1.0) |
| 66 |
self.__reinitPDE.setSolverMethod(solver=LinearPDE.LUMPING) |
| 67 |
self.__h = inf(domain.getSize()) |
| 68 |
self.__smooth = smooth |
| 69 |
self.__n_step=0 |
| 70 |
|
| 71 |
def __advect(self, velocity, dt): |
| 72 |
""" |
| 73 |
Advects the level set function in the presence of a velocity field. |
| 74 |
|
| 75 |
This implementation uses the 2-step Taylor-Galerkin method. |
| 76 |
|
| 77 |
@param velocity: velocity field |
| 78 |
@param dt: time increment |
| 79 |
@return: the advected level set function |
| 80 |
""" |
| 81 |
Y = self.__phi-dt/2.*inner(velocity,grad(self.__phi)) |
| 82 |
self.__PDE.setValue(Y=Y) |
| 83 |
phi_half = self.__PDE.getSolution() |
| 84 |
Y = self.__phi-dt*inner(velocity,grad(phi_half)) |
| 85 |
self.__PDE.setValue(Y=Y) |
| 86 |
phi = self.__PDE.getSolution() |
| 87 |
print "LevelSet: Advection step done" |
| 88 |
return phi |
| 89 |
|
| 90 |
def __reinitialise(self): |
| 91 |
""" |
| 92 |
Reinitializes the level set. |
| 93 |
|
| 94 |
It solves the PDE... |
| 95 |
|
| 96 |
@return: reinitialized level set |
| 97 |
""" |
| 98 |
phi=self.__phi |
| 99 |
s = sign(phi.interpolate(Function(self.__domain))) |
| 100 |
g=grad(phi) |
| 101 |
w = s*g/(length(g)+1e-10) |
| 102 |
dtau = 0.5*self.__h |
| 103 |
iter =0 |
| 104 |
mask = whereNegative(abs(phi)-1.*self.__h) |
| 105 |
self.__reinitPDE.setValue(q=mask, r=phi) |
| 106 |
g=grad(phi) |
| 107 |
while (iter<=self.__reinit_max): |
| 108 |
Y = phi+(dtau/2.)*(s-inner(w,g)) |
| 109 |
self.__reinitPDE.setValue(Y = Y) |
| 110 |
phi_half = self.__reinitPDE.getSolution() |
| 111 |
Y = phi+dtau*(s-inner(w,grad(phi_half))) |
| 112 |
self.__reinitPDE.setValue(Y = Y) |
| 113 |
phi = self.__reinitPDE.getSolution() |
| 114 |
g=grad(phi) |
| 115 |
error = Lsup(length(g)*whereNegative(abs(phi.interpolate(g.getFunctionSpace()))-3.0*self.__h)) |
| 116 |
print "LevelSet:reinitialization: iteration :", iter, " error:", error |
| 117 |
iter +=1 |
| 118 |
return phi |
| 119 |
|
| 120 |
def getTimeStepSize(self,velocity): |
| 121 |
""" |
| 122 |
Returns a new C{dt} for a given C{velocity} using the Courant condition. |
| 123 |
|
| 124 |
@param velocity: velocity field |
| 125 |
""" |
| 126 |
self.__velocity=velocity |
| 127 |
dt=0.5*self.__h/sup(length(velocity)) |
| 128 |
return dt |
| 129 |
|
| 130 |
def update(self,dt): |
| 131 |
""" |
| 132 |
Sets a new velocity and updates the level set function. |
| 133 |
|
| 134 |
@param dt: time step forward |
| 135 |
""" |
| 136 |
self.__phi=self.__advect(self.__velocity, dt) |
| 137 |
self.__n_step+=1 |
| 138 |
if self.__n_step%self.__reinit_each ==0: self.__phi = self.__reinitialise() |
| 139 |
return self.__phi |
| 140 |
|
| 141 |
def update_phi(self, velocity, dt): |
| 142 |
""" |
| 143 |
Updates C{phi} under the presence of a velocity field. |
| 144 |
|
| 145 |
If dt is small this call is equivalent to:: |
| 146 |
|
| 147 |
dt=LevelSet.getTimeStepSize(velocity) |
| 148 |
phi=LevelSet.update(dt) |
| 149 |
|
| 150 |
otherwise substepping is used. |
| 151 |
|
| 152 |
@param velocity: velocity field |
| 153 |
@param dt: time step forward |
| 154 |
""" |
| 155 |
dt2=self.getTimeStepSize(velocity) |
| 156 |
n=math.ceil(dt/dt2) |
| 157 |
dt_new=dt/n |
| 158 |
for i in range(n): |
| 159 |
phi=self.update(dt_new) |
| 160 |
t+=dt_new |
| 161 |
return phi |
| 162 |
|
| 163 |
|
| 164 |
def getVolume(self): |
| 165 |
""" |
| 166 |
Returns the volume of the M{phi(x)<0} region. |
| 167 |
""" |
| 168 |
return integrate(whereNegative(self.__phi.interpolate(Function(self.__domain)))) |
| 169 |
|
| 170 |
def getSurface(self,rel_width_factor=0.5): |
| 171 |
""" |
| 172 |
Returns a mask for the M{phi(x)=1} region |
| 173 |
|
| 174 |
@param rel_width_factor: relative width of region around zero contour |
| 175 |
""" |
| 176 |
return whereNegative(abs(self.__phi)-rel_width_factor*self.__h) |
| 177 |
|
| 178 |
def getH(self): |
| 179 |
""" |
| 180 |
Returns the mesh size. |
| 181 |
""" |
| 182 |
return self.__h |
| 183 |
|
| 184 |
def getDomain(self): |
| 185 |
""" |
| 186 |
Returns the domain. |
| 187 |
""" |
| 188 |
return self.__domain |
| 189 |
|
| 190 |
def getLevelSetFunction(self): |
| 191 |
""" |
| 192 |
Returns the level set function. |
| 193 |
""" |
| 194 |
return self.__phi |
| 195 |
|
| 196 |
def update_parameter_sharp(self, param_neg=-1, param_pos=1, phi=None): |
| 197 |
""" |
| 198 |
Creates a function with C{param_neg} where C{phi<0} and C{param_pos} |
| 199 |
where C{phi>0} (no smoothing). |
| 200 |
|
| 201 |
@param param_neg: value of parameter on the negative side (phi<0) |
| 202 |
@param param_pos: value of parameter on the positive side (phi>0) |
| 203 |
@param phi: level set function to be used. If not present the current |
| 204 |
level set is used. |
| 205 |
""" |
| 206 |
mask_neg = whereNegative(self.__phi) |
| 207 |
mask_pos = whereNonNegative(self.__phi) |
| 208 |
param = param_pos*mask_pos + param_neg*mask_neg |
| 209 |
return param |
| 210 |
|
| 211 |
def update_parameter(self, param_neg=-1, param_pos=1, phi=None, smoothing_width=None): |
| 212 |
""" |
| 213 |
Creates a smoothed function with C{param_neg} where C{phi<0} and |
| 214 |
C{param_pos} where C{phi>0} which is smoothed over a length |
| 215 |
C{smoothing_width} across the interface. |
| 216 |
|
| 217 |
@param smoothing_width: width of the smoothing zone relative to mesh size. |
| 218 |
If not present the initial value of C{smooth} is |
| 219 |
used. |
| 220 |
""" |
| 221 |
if smoothing_width==None: smoothing_width = self.__smooth |
| 222 |
if phi==None: phi=self.__phi |
| 223 |
s=self.__makeInterface(phi,smoothing_width) |
| 224 |
return ((param_pos-param_neg)*s+param_pos+param_neg)/2 |
| 225 |
|
| 226 |
def __makeInterface(self,phi,smoothing_width): |
| 227 |
""" |
| 228 |
Creates a smooth interface from -1 to 1 over the length |
| 229 |
M{2*h*smoothing_width} where -1 is used where the level set is negative |
| 230 |
and 1 where the level set is 1. |
| 231 |
""" |
| 232 |
s=smoothing_width*self.__h |
| 233 |
phi_on_h=interpolate(phi,Function(self.__domain)) |
| 234 |
mask_neg = whereNonNegative(-s-phi_on_h) |
| 235 |
mask_pos = whereNonNegative(phi_on_h-s) |
| 236 |
mask_interface = 1.-mask_neg-mask_pos |
| 237 |
interface=phi_on_h/s |
| 238 |
return - mask_neg + mask_pos + mask_interface * interface |
| 239 |
|
| 240 |
def makeCharacteristicFunction(self, contour=0, phi=None, positiveSide=True, smoothing_width=None): |
| 241 |
""" |
| 242 |
Makes a smooth characteristic function of the region C{phi(x)>contour} if |
| 243 |
C{positiveSide} and C{phi(x)<contour} otherwise. |
| 244 |
|
| 245 |
@param phi: level set function to be used. If not present the current |
| 246 |
level set is used. |
| 247 |
@param smoothing_width: width of the smoothing zone relative to mesh size. |
| 248 |
If not present the initial value of C{smooth} is |
| 249 |
used. |
| 250 |
""" |
| 251 |
if phi==None: phi=self.__phi |
| 252 |
if smoothing_width == None: smoothing_width=self.__smooth |
| 253 |
s=self.__makeInterface(phi=phi-contour,smoothing_width=smoothing_width) |
| 254 |
if positiveSide: |
| 255 |
return (1+s)/2 |
| 256 |
else: |
| 257 |
return (1-s)/2 |
| 258 |
|
| 259 |
def setTolerance(self,tolerance=1e-3): |
| 260 |
self.__PDE.setTolerance(tolerance) |
| 261 |
self.__reinitPDE.setTolerance(tolerance) |
| 262 |
|
| 263 |
|
| 264 |
class LevelSet2(object): |
| 265 |
def __init__(self,phi,reinit_max=10,reinit_each=2,smooth=2.): |
| 266 |
""" |
| 267 |
Initializes the model. |
| 268 |
""" |
| 269 |
self.__domain = phi.getDomain() |
| 270 |
x=self.__domain.getX() |
| 271 |
diam=0 |
| 272 |
for i in range(self.__domain.getDim()): |
| 273 |
xi=x[i] |
| 274 |
diam+=(inf(xi)-sup(xi))**2 |
| 275 |
self.__diam=sqrt(diam) |
| 276 |
self.__h = sup(Function(self.__domain).getSize()) |
| 277 |
self.__reinit_each=reinit_each |
| 278 |
self.__reinit_max = reinit_max |
| 279 |
self.__smooth = smooth |
| 280 |
self.__phi = phi |
| 281 |
self.__update_count=0 |
| 282 |
self.velocity = None |
| 283 |
|
| 284 |
#================================================== |
| 285 |
self.__FC=True |
| 286 |
if self.__FC: |
| 287 |
self.__fc=TransportPDE(self.__domain,num_equations=1,useBackwardEuler=False) |
| 288 |
# self.__fc.setReducedOrderOn() |
| 289 |
self.__fc.setValue(M=Scalar(1.,Function(self.__domain))) |
| 290 |
self.__fc.setInitialSolution(phi+self.__diam) |
| 291 |
else: |
| 292 |
self.__fcpde = LinearPDE(self.__domain,numEquations=1, numSolutions=1) |
| 293 |
self.__fcpde.setSolverMethod(solver=LinearPDE.LUMPING) |
| 294 |
self.__fcpde.setReducedOrderOn() |
| 295 |
self.__fcpde.setSymmetryOn() |
| 296 |
self.__fcpde.setValue(D=1.) |
| 297 |
|
| 298 |
#======================================= |
| 299 |
self.__reinitFC=False |
| 300 |
if self.__reinitFC: |
| 301 |
self.__reinitfc=TransportPDE(self.__domain,num_equations=1,useBackwardEuler=True) |
| 302 |
self.__reinitfc.setValue(M=Scalar(1.,Function(self.__domain))) |
| 303 |
self.__reinitfc.setTolerance(1.e-5) |
| 304 |
else: |
| 305 |
self.__reinitPde = LinearPDE(self.__domain,numEquations=1, numSolutions=1) |
| 306 |
# self.__reinitPde.setSolverMethod(solver=LinearPDE.LUMPING) |
| 307 |
self.__reinitPde.setReducedOrderOn() |
| 308 |
self.__reinitPde.setSymmetryOn() |
| 309 |
self.__reinitPde.setValue(D=1.) |
| 310 |
self.__reinitPde.setTolerance(1.e-2) |
| 311 |
# self.__reinitPde.setSolverMethod(preconditioner=LinearPDE.ILU0) |
| 312 |
#======================================= |
| 313 |
self.__updateInterface() |
| 314 |
print "phi range:",inf(phi), sup(phi) |
| 315 |
|
| 316 |
def setFixedRegion(self,mask,contour=0): |
| 317 |
q=whereNonPositive(abs(self.__phi-contour)-self.__h)*mask |
| 318 |
# self.__fc.setValue(q=q) |
| 319 |
self.__reinitPde.setValue(q=q) |
| 320 |
|
| 321 |
|
| 322 |
def __updateInterface(self): |
| 323 |
self.__smoothed_char=self.__makeInterface(self.__phi) |
| 324 |
|
| 325 |
|
| 326 |
def update(self,dt): |
| 327 |
""" |
| 328 |
Sets a new velocity and updates the level set function. |
| 329 |
|
| 330 |
@param dt: time step forward |
| 331 |
""" |
| 332 |
if self.__FC: |
| 333 |
self.__phi=self.__fc.solve(dt, verbose=False)-self.__diam |
| 334 |
else: |
| 335 |
self.__fcpde.setValue(Y = self.__phi-(dt/2.)*inner(self.velocity,grad(self.__phi))) |
| 336 |
phi_half = self.__fcpde.getSolution() |
| 337 |
self.__fcpde.setValue(Y = self.__phi-dt*inner(self.velocity,grad(phi_half))) |
| 338 |
self.__phi= self.__fcpde.getSolution() |
| 339 |
self.__update_count += 1 |
| 340 |
if self.__update_count%self.__reinit_each == 0: |
| 341 |
self.__phi=self.__reinitialise(self.__phi) |
| 342 |
if self.__FC: |
| 343 |
self.__fc.setInitialSolution(self.__phi+self.__diam) |
| 344 |
self.__updateInterface() |
| 345 |
|
| 346 |
def setTolerance(self,tolerance=1e-3): |
| 347 |
if self.__FC: |
| 348 |
self.__fc.setTolerance(tolerance) |
| 349 |
|
| 350 |
def __reinitialise(self,phi): |
| 351 |
print "reintialization started:" |
| 352 |
s=self.__makeInterface(phi,1.) |
| 353 |
print "phi range:",inf(phi), sup(phi) |
| 354 |
g=grad(phi) |
| 355 |
w=s*g/(length(g)+EPSILON) |
| 356 |
#======================================================= |
| 357 |
# # positive part: |
| 358 |
# phi_p=wherePositive(phi)*phi |
| 359 |
# self.__reinitfc.setInitialSolution(phi_p) |
| 360 |
# self.__reinitfc.setValue(C=-wherePositive(s)*w,Y=wherePositive(s)*s,q=whereNonPositive(phi)) |
| 361 |
# dtau=self.__h |
| 362 |
# print "step size: dt (pos)= ",dtau |
| 363 |
# print "phi_p range:",inf(phi_p), sup(phi_p) |
| 364 |
# iter=0 |
| 365 |
# while (iter<=self.__reinit_max): |
| 366 |
# phi_p=self.__reinitfc.solve(dtau) |
| 367 |
# print "phi_p range:",inf(phi_p), sup(phi_p) |
| 368 |
# iter+=1 |
| 369 |
# # negative part: |
| 370 |
# phi_n=-whereNegative(phi)*phi |
| 371 |
# self.__reinitfc.setInitialSolution(phi_n) |
| 372 |
# self.__reinitfc.setValue(C=-whereNegative(s)*w,Y=-whereNegative(s)*s,q=whereNonNegative(phi)) |
| 373 |
# # dtau=self.__reinitfc.getSafeTimeStepSize() |
| 374 |
# dtau=self.__h |
| 375 |
# print "step size: dt (neg)= ",dtau |
| 376 |
# print "phi_n range:",inf(phi_n), sup(phi_n) |
| 377 |
# iter=0 |
| 378 |
# while (iter<=self.__reinit_max): |
| 379 |
# phi_n=self.__reinitfc.solve(dtau) |
| 380 |
# print "phi_n range:",inf(phi_n), sup(phi_n) |
| 381 |
# iter+=1 |
| 382 |
# phi=phi_p-phi_n |
| 383 |
# print "phi range:",inf(phi), sup(phi) |
| 384 |
# print "reintialization completed." |
| 385 |
# return phi |
| 386 |
|
| 387 |
#======================================================= |
| 388 |
if self.__reinitFC: |
| 389 |
self.__reinitfc.setValue(C=-w,Y=s) |
| 390 |
self.__reinitfc.setInitialSolution(phi+self.__diam) |
| 391 |
dtau=self.__reinitfc.getSafeTimeStepSize() |
| 392 |
# dtau = 0.5*inf(Function(self.__domain).getSize()) |
| 393 |
else: |
| 394 |
dtau = 0.5*inf(Function(self.__domain).getSize()) |
| 395 |
print "step size: dt = ",dtau,inf(abs(phi.interpolate(Function(self.__domain)))/abs(s-inner(w,grad(phi)))) |
| 396 |
iter =0 |
| 397 |
# self.__reinitPde.setValue(q=whereNegative(abs(phi)-2*self.__h), r=phi) |
| 398 |
# self.__reinitPde.setValue(r=phi) |
| 399 |
while (iter<=self.__reinit_max): |
| 400 |
phi_old=phi |
| 401 |
if self.__reinitFC: |
| 402 |
phi = self.__reinitfc.solve(dtau)-self.__diam |
| 403 |
else: |
| 404 |
self.__reinitPde.setValue(Y = phi+(dtau/2.)*(s-inner(w,grad(phi)))) |
| 405 |
phi_half = self.__reinitPde.getSolution() |
| 406 |
self.__reinitPde.setValue(Y = phi+dtau*(s-inner(w,grad(phi_half)))) |
| 407 |
# g=grad(phi) |
| 408 |
# S=inner(w,grad(phi)) |
| 409 |
# self.__reinitPde.setValue(Y = phi+dtau*(s-S),X=dtau**2/2*w*(s-S)) |
| 410 |
phi = self.__reinitPde.getSolution() |
| 411 |
change = Lsup(phi-phi_old)/self.__diam |
| 412 |
print "phi range:",inf(phi), sup(phi) |
| 413 |
print "iteration :", iter, " change:", change |
| 414 |
iter +=1 |
| 415 |
print "reintialization completed." |
| 416 |
return phi |
| 417 |
|
| 418 |
def createParameter(self,value_negative=-1.,value_positive=1): |
| 419 |
out = (value_negative+value_positive)/2. + self.__smoothed_char * ((value_positive-value_negative)/2.) |
| 420 |
return out |
| 421 |
|
| 422 |
|
| 423 |
#=========== things from here onwards are not used nor tested: =========== |
| 424 |
|
| 425 |
|
| 426 |
def getCharacteristicFunction(self): |
| 427 |
return self.__smoothed_char |
| 428 |
|
| 429 |
|
| 430 |
def RK2(self,L): |
| 431 |
k0=L(phi) |
| 432 |
phi_1=phi+dt*k0 |
| 433 |
k1=L(phi_1) |
| 434 |
phi=phi+dt/2*(k0+k1) |
| 435 |
|
| 436 |
def RK3(self,L): |
| 437 |
k0=L(phi) |
| 438 |
phi_1=phi+dt*L(phi) |
| 439 |
k1=L(phi_1) |
| 440 |
phi_2=phi+dt/4*(k0+k1) |
| 441 |
k2=L(phi_2) |
| 442 |
phi=phi+dt/6*(k0+4*k2+K1) |
| 443 |
|
| 444 |
def TG(self,L): |
| 445 |
k0=L(phi) |
| 446 |
phi_1=phi+dt/2*k0 |
| 447 |
k1=L(phi_1) |
| 448 |
phi=phi+dt*k1 |
| 449 |
|
| 450 |
|
| 451 |
def __reinitialise_old(self): |
| 452 |
#============================================= |
| 453 |
f=0.1 |
| 454 |
f2=0.3 |
| 455 |
h=self.__h |
| 456 |
fs=h.getFunctionSpace() |
| 457 |
vol=self.getVolumeOfNegativeDomain() |
| 458 |
self.__reinitPde.setValue(D=1.0) |
| 459 |
#============ |
| 460 |
grad_phi=grad(self.__phi,fs) |
| 461 |
len_grad_phi=length(grad_phi) |
| 462 |
w=grad_phi/len_grad_phi |
| 463 |
|
| 464 |
# s=wherePositive(self. __makeInterface(2.))-whereNegative(self. __makeInterface(2.)) |
| 465 |
s=self. __makeInterface(2.) |
| 466 |
# s2=whereNegative(len_grad_phi-f2)*(1-(1-len_grad_phi/f2)**2)+(1-whereNegative(len_grad_phi-f2)) |
| 467 |
s2=1. |
| 468 |
s=s*s2 |
| 469 |
dtau=f*inf(h/abs(s)) |
| 470 |
|
| 471 |
#======================== |
| 472 |
diff=1. |
| 473 |
d=1. |
| 474 |
c =0 |
| 475 |
#============================================== |
| 476 |
TVD=integrate(length(grad_phi)) |
| 477 |
print "initial range ",inf(self.__phi),sup(self.__phi),"error:",Lsup(1.-len_grad_phi),"volume =",vol,TVD |
| 478 |
# saveVTK("test.%s.xml"%c,l=length(grad(self.__phi,fs))-1,s=s,phi=self.__phi) |
| 479 |
|
| 480 |
dtau=f*inf(h/abs(s)) |
| 481 |
while c < self.__reinit_max: # and abs(diff) >= 0.01: |
| 482 |
# |
| 483 |
grad_phi=grad(self.__phi,fs) |
| 484 |
len_grad_phi=length(grad_phi) |
| 485 |
# |
| 486 |
# s=self.__makeInterface(1.) |
| 487 |
# s2=whereNegative(len_grad_phi-f2)*(1-(1-len_grad_phi/f2)**2)+(1-whereNegative(len_grad_phi-f2)) |
| 488 |
# s=s*s2 |
| 489 |
# phi_on_h=interpolate(self.__phi,fs) |
| 490 |
# self.__reinitPde.setValue(Y =self.__phi+dtau/2*s*(1.-len_grad_phi)) |
| 491 |
# phi_half=self.__reinitPde.getSolution() |
| 492 |
# self.__reinitPde.setValue(Y =self.__phi+dtau*s*(1.-length(grad(phi_half,fs)))) |
| 493 |
# self.__reinitPde.setValue(Y =self.__phi, Y_reduced=dtau*s*(1.-inner(w,grad(phi_half,ReducedFunction(self.__domain))))) |
| 494 |
# self.__reinitPde.setValue(Y =self.__phi+dtau*(s-inner(w,grad_phi)),X=dtau/2*h*(s-inner(w,grad_phi))*w) |
| 495 |
# self.__reinitPde.setValue(Y =self.__phi+dtau*s*(1.-len_grad_phi),X=-dtau/2*h*s**2*grad_phi) |
| 496 |
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) |
| 497 |
|
| 498 |
# 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) |
| 499 |
self.__phi, previous = self.__reinitPde.getSolution(), self.__phi |
| 500 |
# self.__updateInterface() |
| 501 |
|
| 502 |
vol,vol_old=self.getVolumeOfNegativeDomain(),vol |
| 503 |
diff=(vol-vol_old) |
| 504 |
r=Lsup(length(grad(self.__phi))-1.) |
| 505 |
TVD=integrate(length(grad(self.__phi,fs))) |
| 506 |
print "iteration :", c, "range ",inf(self.__phi),sup(self.__phi),"error :",r,"volume change:",diff,TVD |
| 507 |
# saveVTK("test.%s.xml"%(c+1),l=length(grad(self.__phi,fs)),s=s,phi=self.__phi,v=grad(self.__phi,fs)) |
| 508 |
c += 1 |
| 509 |
return |
| 510 |
#============================================== |
| 511 |
f2=0.7 |
| 512 |
h=self.__h |
| 513 |
fs=h.getFunctionSpace() |
| 514 |
vol=self.getVolumeOfNegativeDomain() |
| 515 |
s=abs(self. __makeInterface(2.)) |
| 516 |
grad_phi=grad(self.__phi,fs) |
| 517 |
len_grad_phi=length(grad_phi) |
| 518 |
|
| 519 |
#---- |
| 520 |
# aphi_on_h=abs(interpolate(self.__phi,self.__h.getFunctionSpace())) |
| 521 |
# q=2*self.__h |
| 522 |
# mask_neg = whereNegative(aphi_on_h-q) |
| 523 |
# mask_pos = wherePositive(aphi_on_h-q-self.__h) |
| 524 |
# mask_interface = 1.-mask_neg-mask_pos |
| 525 |
# s2=mask_neg + mask_interface * (q+self.__h-aphi_on_h)/self.__h |
| 526 |
#---- |
| 527 |
m=whereNonPositive(len_grad_phi-f2) |
| 528 |
s2=m*len_grad_phi/f2+(1.-m) |
| 529 |
#---- |
| 530 |
# s2=1. |
| 531 |
#---- |
| 532 |
self.__reinitPde.setValue(D=(1-s)/h**2,Y=(1-s)/h**2*self.__phi) |
| 533 |
self.__reinitPde.setValue(A=(s2*len_grad_phi+(1.-s2))*kronecker(3),X=grad_phi) |
| 534 |
# self.__reinitPde.setValue(A=(len_grad_phi-1)/len_grad_phi*kronecker(3)) |
| 535 |
# self.__reinitPde.setValue(A=kronecker(3), X=grad_phi/len_grad_phi) |
| 536 |
self.__phi = self.__reinitPde.getSolution() |
| 537 |
r=Lsup(length(grad(self.__phi))-1.) |
| 538 |
vol,vol_old=self.getVolumeOfNegativeDomain(),vol |
| 539 |
diff=(vol-vol_old)/vol |
| 540 |
print "iteration :", inf(self.__phi),sup(self.__phi),r,diff |
| 541 |
# saveVTK("test.%s.xml"%0,l=length(grad(self.__phi,fs)),s=s,phi=self.__phi,v=grad(self.__phi,fs),s2=s2) |
| 542 |
return |
| 543 |
#============================================= |
| 544 |
|
| 545 |
#============ |
| 546 |
|
| 547 |
|
| 548 |
def getVolumeOfNegativeDomain(self): |
| 549 |
""" |
| 550 |
Returns the current volume of domain with phi<0. |
| 551 |
""" |
| 552 |
return integrate((1.-self.__makeInterface(1.))/2.) |
| 553 |
|
| 554 |
|
| 555 |
def getBoundingBoxOfNegativeDomain(self): |
| 556 |
""" |
| 557 |
Returns the height of the region with phi<0. |
| 558 |
""" |
| 559 |
fs=self.__h.getFunctionSpace() |
| 560 |
mask_phi1=wherePositive(interpolate(self.__phi,fs)) |
| 561 |
mask_phi2=wherePositive(self.__phi) |
| 562 |
x1=fs.getX() |
| 563 |
x2=self.__domain.getX() |
| 564 |
out=[] |
| 565 |
for i in range(fs.getDim()): |
| 566 |
x1_i=x1[i] |
| 567 |
x2_i=x2[i] |
| 568 |
d=2*(sup(x2_i)-inf(x2_i)) |
| 569 |
offset1=d*mask_phi1 |
| 570 |
offset2=d*mask_phi2 |
| 571 |
out.append((min(inf(x1_i+offset1),inf(x2_i+offset2)),max(sup(x1_i-offset1),sup(x2_i-offset2)))) |
| 572 |
return tuple(out) |
| 573 |
|
Parent Directory
| 
Revision Log