# Contents of /trunk/escript/py_src/levelset.py

Revision 3432 - (show annotations)
Fri Jan 7 01:32:07 2011 UTC (10 years, 4 months ago) by jfenwick
File MIME type: text/x-python
File size: 21862 byte(s)
```Made import statements a bit more specific to clean up the epydoc
```
 1 2 ######################################################## 3 # 4 # Copyright (c) 2003-2010 by University of Queensland 5 # Earth Systems Science Computational Center (ESSCC) 6 7 # 8 # Primary Business: Queensland, Australia 9 # Licensed under the Open Software License version 3.0 10 11 # 12 ######################################################## 13 14 __copyright__="""Copyright (c) 2003-2010 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 20 __url__= 21 22 #from esys.escript import * 23 from esys.escript.linearPDEs import LinearPDE, SingleTransportPDE 24 from esys.escript.pdetools import Projector 25 26 import math 27 28 USE_OLD_VERSION=False 29 USE_OLD_VERSION_REINIT=False 30 31 32 class LevelSet: 33 """ 34 The level set method tracking an interface defined by the zero contour of the 35 level set function phi which defines the signed distance of a point x from the 36 interface. The contour phi(x)=0 defines the interface. 37 38 It is assumed that phi(x)<0 defines the volume of interest, 39 phi(x)>0 the outside world. 40 """ 41 def __init__(self,phi,reinit_max=10,reinitialize_after=20,smooth=2., useReducedOrder=False): 42 """ 43 Sets up the level set method. 44 45 :param phi: the initial level set function 46 :param reinit_max: maximum number of reinitialization steps 47 :param reinitialize_after: ``phi`` is reinitialized every ``reinit_after`` step 48 :param smooth: smoothing width 49 """ 50 self.__domain = phi.getDomain() 51 self.__phi = phi 52 53 if USE_OLD_VERSION: 54 self.__PDE = LinearPDE(self.__domain) 55 if useReducedOrder: self.__PDE.setReducedOrderOn() 56 self.__PDE.setSolverMethod(solver=LinearPDE.PCG) 57 self.__PDE.setValue(D=1.0) 58 else: 59 self.__PDE=SingleTransportPDE(self.__domain) 60 if useReducedOrder: self.__PDE.setReducedOrderOn() 61 self.__PDE.setValue(M=1.0) 62 63 if USE_OLD_VERSION_REINIT: 64 self.__reinitPDE = LinearPDE(self.__domain, numEquations=1) 65 self.__reinitPDE.getSolverOption().setSolverMethod(solver=LinearPDE.LUMPING) 66 if useReducedOrder: self.__reinitPDE.setReducedOrderOn() 67 self.__reinitPDE.setValue(D=1.0) 68 else: 69 self.__reinitPDE=SingleTransportPDE(self.__domain) 70 if useReducedOrder: self.__reinitPDE.setReducedOrderOn() 71 self.__reinitPDE.setValue(M=1.0) 72 73 # revise: 74 self.__reinit_max = reinit_max 75 self.__reinit_after = reinitialize_after 76 self.__h = inf(self.__domain.getSize()) 77 self.__smooth = smooth 78 self.__n_step=0 79 80 def getH(self): 81 """ 82 Returns the mesh size. 83 """ 84 return self.__h 85 86 def getDomain(self): 87 """ 88 Returns the domain. 89 """ 90 return self.__domain 91 92 def getAdvectionSolverOptions(self): 93 """ 94 Returns the solver options for the interface advective. 95 """ 96 return self.__PDE.getSolverOptions() 97 98 def getReinitializationSolverOptions(self): 99 """ 100 Returns the options of the solver for the reinitialization 101 """ 102 return self.__reinitPDE.getSolverOption() 103 104 def getLevelSetFunction(self): 105 """ 106 Returns the level set function. 107 """ 108 return self.__phi 109 110 def getTimeStepSize(self,velocity): 111 """ 112 Returns a new ``dt`` for a given ``velocity`` using the Courant condition. 113 114 :param velocity: velocity field 115 """ 116 if USE_OLD_VERSION: 117 self.__velocity=velocity 118 dt=0.5*self.__h/sup(length(velocity)) 119 else: 120 self.__PDE.setValue(C=-velocity) 121 dt=self.__PDE.getSafeTimeStepSize() 122 123 return dt 124 125 def update(self,dt): 126 """ 127 Sets a new velocity and updates the level set function. 128 129 :param dt: time step forward 130 """ 131 self.__phi=self.__advect(dt) 132 self.__n_step+=1 133 if self.__n_step%self.__reinit_after ==0: self.__phi = self.__reinitialise() 134 return self.__phi 135 136 def update_phi(self, velocity, dt): 137 """ 138 Updates ``phi`` under the presence of a velocity field. 139 140 If dt is small this call is equivalent to:: 141 142 dt=LevelSet.getTimeStepSize(velocity) 143 phi=LevelSet.update(dt) 144 145 otherwise substepping is used. 146 147 :param velocity: velocity field 148 :param dt: time step forward 149 """ 150 dt2=self.getTimeStepSize(velocity) 151 n=int(math.ceil(dt/dt2)) 152 dt_new=dt/n 153 for i in range(n): 154 phi=self.update(dt_new) 155 return phi 156 157 def __advect(self, dt): 158 """ 159 Advects the level set function in the presence of a velocity field. 160 161 :param dt: time increment 162 :return: the advected level set function 163 """ 164 if USE_OLD_VERSION: 165 velocity=self.__velocity 166 Y = self.__phi-dt/2.*inner(velocity,grad(self.__phi)) 167 self.__PDE.setValue(Y=Y) 168 phi_half = self.__PDE.getSolution() 169 Y = self.__phi-dt*inner(velocity,grad(phi_half)) 170 self.__PDE.setValue(Y=Y) 171 phi = self.__PDE.getSolution() 172 else: 173 phi = self.__PDE.getSolution(dt, self.__phi) 174 print "LevelSet: Advection step done" 175 return phi 176 177 #============================================================================================== 178 def __reinitialise(self): 179 """ 180 Reinitializes the level set. 181 182 It solves the PDE... 183 184 :return: reinitialized level set 185 """ 186 phi=self.__phi 187 s = sign(phi.interpolate(Function(self.__domain))) 188 g=grad(phi) 189 w = s*g/(length(g)+1e-10) 190 dtau = 0.5*self.__h 191 iter =0 192 mask = whereNegative(abs(phi)-1.*self.__h) 193 self.__reinitPDE.setValue(q=mask, r=phi) 194 g=grad(phi) 195 while (iter<=self.__reinit_max): 196 Y = phi+(dtau/2.)*(s-inner(w,g)) 197 self.__reinitPDE.setValue(Y = Y) 198 phi_half = self.__reinitPDE.getSolution() 199 Y = phi+dtau*(s-inner(w,grad(phi_half))) 200 self.__reinitPDE.setValue(Y = Y) 201 phi = self.__reinitPDE.getSolution() 202 g=grad(phi) 203 error = Lsup(length(g)*whereNegative(abs(phi.interpolate(g.getFunctionSpace()))-3.0*self.__h)) 204 print "LevelSet:reinitialization: iteration :", iter, " error:", error 205 iter +=1 206 return phi 207 208 209 def getVolume(self): 210 """ 211 Returns the volume of the *phi(x)<0* region. 212 """ 213 return integrate(whereNegative(self.__phi.interpolate(Function(self.__domain)))) 214 215 216 def update_parameter_sharp(self, param_neg=-1, param_pos=1, phi=None): 217 """ 218 Creates a function with ``param_neg`` where ``phi<0`` and ``param_pos`` 219 where ``phi>0`` (no smoothing). 220 221 :param param_neg: value of parameter on the negative side (phi<0) 222 :param param_pos: value of parameter on the positive side (phi>0) 223 :param phi: level set function to be used. If not present the current 224 level set is used. 225 """ 226 mask_neg = whereNegative(self.__phi) 227 mask_pos = whereNonNegative(self.__phi) 228 param = param_pos*mask_pos + param_neg*mask_neg 229 return param 230 231 def update_parameter(self, param_neg=-1, param_pos=1, phi=None, smoothing_width=None): 232 """ 233 Creates a smoothed function with ``param_neg`` where ``phi<0`` and 234 ``param_pos`` where ``phi>0`` which is smoothed over a length 235 ``smoothing_width`` across the interface. 236 237 :param smoothing_width: width of the smoothing zone relative to mesh size. 238 If not present the initial value of ``smooth`` is 239 used. 240 """ 241 if smoothing_width==None: smoothing_width = self.__smooth 242 if phi==None: phi=self.__phi 243 s=self.__makeInterface(phi,smoothing_width) 244 return ((param_pos-param_neg)*s+param_pos+param_neg)/2 245 246 def __makeInterface(self,phi,smoothing_width): 247 """ 248 Creates a smooth interface from -1 to 1 over the length 249 *2*h*smoothing_width* where -1 is used where the level set is negative 250 and 1 where the level set is 1. 251 """ 252 s=smoothing_width*self.__h 253 phi_on_h=interpolate(phi,Function(self.__domain)) 254 mask_neg = whereNonNegative(-s-phi_on_h) 255 mask_pos = whereNonNegative(phi_on_h-s) 256 mask_interface = 1.-mask_neg-mask_pos 257 interface=phi_on_h/s 258 return - mask_neg + mask_pos + mask_interface * interface 259 260 def makeCharacteristicFunction(self, contour=0, phi=None, positiveSide=True, smoothing_width=None): 261 """ 262 Makes a smooth characteristic function of the region ``phi(x)>contour`` if 263 ``positiveSide`` and ``phi(x)= 0.01: 498 # 499 grad_phi=grad(self.__phi,fs) 500 len_grad_phi=length(grad_phi) 501 # 502 # s=self.__makeInterface(1.) 503 # s2=whereNegative(len_grad_phi-f2)*(1-(1-len_grad_phi/f2)**2)+(1-whereNegative(len_grad_phi-f2)) 504 # s=s*s2 505 # phi_on_h=interpolate(self.__phi,fs) 506 # self.__reinitPde.setValue(Y =self.__phi+dtau/2*s*(1.-len_grad_phi)) 507 # phi_half=self.__reinitPde.getSolution() 508 # self.__reinitPde.setValue(Y =self.__phi+dtau*s*(1.-length(grad(phi_half,fs)))) 509 # self.__reinitPde.setValue(Y =self.__phi, Y_reduced=dtau*s*(1.-inner(w,grad(phi_half,ReducedFunction(self.__domain))))) 510 # self.__reinitPde.setValue(Y =self.__phi+dtau*(s-inner(w,grad_phi)),X=dtau/2*h*(s-inner(w,grad_phi))*w) 511 # self.__reinitPde.setValue(Y =self.__phi+dtau*s*(1.-len_grad_phi),X=-dtau/2*h*s**2*grad_phi) 512 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) 513 514 # 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) 515 self.__phi, previous = self.__reinitPde.getSolution(), self.__phi 516 # self.__updateInterface() 517 518 vol,vol_old=self.getVolumeOfNegativeDomain(),vol 519 diff=(vol-vol_old) 520 r=Lsup(length(grad(self.__phi))-1.) 521 TVD=integrate(length(grad(self.__phi,fs))) 522 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)) 524 c += 1 525 return 526 #============================================== 527 f2=0.7 528 h=self.__h 529 fs=h.getFunctionSpace() 530 vol=self.getVolumeOfNegativeDomain() 531 s=abs(self. __makeInterface(2.)) 532 grad_phi=grad(self.__phi,fs) 533 len_grad_phi=length(grad_phi) 534 535 #---- 536 # aphi_on_h=abs(interpolate(self.__phi,self.__h.getFunctionSpace())) 537 # q=2*self.__h 538 # mask_neg = whereNegative(aphi_on_h-q) 539 # mask_pos = wherePositive(aphi_on_h-q-self.__h) 540 # mask_interface = 1.-mask_neg-mask_pos 541 # s2=mask_neg + mask_interface * (q+self.__h-aphi_on_h)/self.__h 542 #---- 543 m=whereNonPositive(len_grad_phi-f2) 544 s2=m*len_grad_phi/f2+(1.-m) 545 #---- 546 # s2=1. 547 #---- 548 self.__reinitPde.setValue(D=(1-s)/h**2,Y=(1-s)/h**2*self.__phi) 549 self.__reinitPde.setValue(A=(s2*len_grad_phi+(1.-s2))*kronecker(3),X=grad_phi) 550 # self.__reinitPde.setValue(A=(len_grad_phi-1)/len_grad_phi*kronecker(3)) 551 # self.__reinitPde.setValue(A=kronecker(3), X=grad_phi/len_grad_phi) 552 self.__phi = self.__reinitPde.getSolution() 553 r=Lsup(length(grad(self.__phi))-1.) 554 vol,vol_old=self.getVolumeOfNegativeDomain(),vol 555 diff=(vol-vol_old)/vol 556 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) 558 return 559 #============================================= 560 561 #============ 562 563 564 def getVolumeOfNegativeDomain(self): 565 """ 566 Returns the current volume of domain with phi<0. 567 """ 568 return integrate((1.-self.__makeInterface(1.))/2.) 569 570 571 def getBoundingBoxOfNegativeDomain(self): 572 """ 573 Returns the height of the region with phi<0. 574 """ 575 fs=self.__h.getFunctionSpace() 576 mask_phi1=wherePositive(interpolate(self.__phi,fs)) 577 mask_phi2=wherePositive(self.__phi) 578 x1=fs.getX() 579 x2=self.__domain.getX() 580 out=[] 581 for i in range(fs.getDim()): 582 x1_i=x1[i] 583 x2_i=x2[i] 584 d=2*(sup(x2_i)-inf(x2_i)) 585 offset1=d*mask_phi1 586 offset2=d*mask_phi2 587 out.append((min(inf(x1_i+offset1),inf(x2_i+offset2)),max(sup(x1_i-offset1),sup(x2_i-offset2)))) 588 return tuple(out) 589

 ViewVC Help Powered by ViewVC 1.1.26