| 1 |
from esys.escript import * |
| 2 |
from esys.escript.linearPDEs import LinearPDE, TransportPDE |
| 3 |
from esys.escript.pdetools import Projector |
| 4 |
from esys import finley |
| 5 |
import math |
| 6 |
|
| 7 |
|
| 8 |
|
| 9 |
class LevelSet(object): |
| 10 |
def __init__(self,phi,reinitialization_steps_max=10,reinitialization_each=2,relative_smoothing_width=2.): |
| 11 |
""" |
| 12 |
initialize model |
| 13 |
""" |
| 14 |
self.__domain = phi.getDomain() |
| 15 |
x=self.__domain.getX() |
| 16 |
diam=0 |
| 17 |
for i in range(self.__domain.getDim()): |
| 18 |
xi=x[i] |
| 19 |
diam+=(inf(xi)-sup(xi))**2 |
| 20 |
self.__diam=sqrt(diam) |
| 21 |
self.__h = sup(Function(self.__domain).getSize()) |
| 22 |
self.__reinitialization_each=reinitialization_each |
| 23 |
self.__reinitialization_steps_max = reinitialization_steps_max |
| 24 |
self.__relative_smoothing_width = relative_smoothing_width |
| 25 |
self.__phi = phi |
| 26 |
self.__update_count=0 |
| 27 |
self.velocity = None |
| 28 |
|
| 29 |
#================================================== |
| 30 |
self.__fc=TransportPDE(self.__domain,num_equations=1,theta=0.5) |
| 31 |
# self.__fc.setReducedOrderOn() |
| 32 |
self.__fc.setValue(M=Scalar(1.,Function(self.__domain))) |
| 33 |
self.__fc.setInitialSolution(phi+self.__diam) |
| 34 |
|
| 35 |
#======================================= |
| 36 |
self.__reinitPde = LinearPDE(self.__domain,numEquations=1, numSolutions=1) |
| 37 |
self.__reinitPde.setSolverMethod(solver=LinearPDE.LUMPING) |
| 38 |
self.__reinitPde.setReducedOrderOn() |
| 39 |
self.__reinitPde.setSymmetryOn() |
| 40 |
self.__reinitPde.setValue(D=1.) |
| 41 |
# self.__reinitPde.setTolerance(1.e-8) |
| 42 |
# self.__reinitPde.setSolverMethod(preconditioner=LinearPDE.ILU0) |
| 43 |
#======================================= |
| 44 |
self.__updateInterface() |
| 45 |
print "phi range:",inf(phi), sup(phi) |
| 46 |
|
| 47 |
def setFixedRegion(self,mask,contour=0): |
| 48 |
q=whereNonPositive(abs(self.__phi-contour)-self.__h)*mask |
| 49 |
# self.__fc.setValue(q=q) |
| 50 |
self.__reinitPde.setValue(q=q) |
| 51 |
|
| 52 |
def getDomain(self): |
| 53 |
return self.__domain |
| 54 |
|
| 55 |
def getTimeStepSize(self,velocity): |
| 56 |
""" |
| 57 |
returns a new dt for a given velocity using the courant coundition |
| 58 |
""" |
| 59 |
self.velocity=velocity |
| 60 |
self.__fc.setValue(C=-interpolate(velocity,Function(self.__domain))) |
| 61 |
dt=self.__fc.getSafeTimeStepSize() |
| 62 |
return dt |
| 63 |
|
| 64 |
def getLevelSetFunction(self): |
| 65 |
""" |
| 66 |
returns the level set function |
| 67 |
""" |
| 68 |
return self.__phi |
| 69 |
|
| 70 |
def __updateInterface(self): |
| 71 |
self.__smoothed_char=self.__makeInterface(self.__phi) |
| 72 |
|
| 73 |
def __makeInterface(self,phi,smoothing_width=1.): |
| 74 |
""" |
| 75 |
creates a very smooth interface from -1 to 1 over the length 2*h*smoothing_width where -1 is used where the level set is negative |
| 76 |
and 1 where the level set is 1 |
| 77 |
""" |
| 78 |
s=smoothing_width*self.__h |
| 79 |
phi_on_h=interpolate(phi,Function(self.__domain)) |
| 80 |
mask_neg = whereNonNegative(-s-phi_on_h) |
| 81 |
mask_pos = whereNonNegative(phi_on_h-s) |
| 82 |
mask_interface = 1.-mask_neg-mask_pos |
| 83 |
# interface=1.-(phi_on_h-s)**2/(2.*s**3)*(phi_on_h+2*s) # function f with f(s)=1, f'(s)=0, f(-s)=-1, f'(-s)=0, f(0)=0, f'(0)= |
| 84 |
interface=phi_on_h/s |
| 85 |
return - mask_neg + mask_pos + mask_interface * interface |
| 86 |
|
| 87 |
def makeCharacteristicFunction(self, contour=0, positiveSide=True, smoothing_width=1.): |
| 88 |
return self.makeCharacteristicFunctionFromExternalLevelSetFunction(self.__phi,contour,positiveSide,smoothing_width) |
| 89 |
|
| 90 |
def makeCharacteristicFunctionFromExternalLevelSetFunction(self,phi,contour=0, positiveSide=True, smoothing_width=1.): |
| 91 |
s=self.__makeInterface(phi-contour,smoothing_width) |
| 92 |
if positiveSide: |
| 93 |
return (1+s)/2 |
| 94 |
else: |
| 95 |
return (1-s)/2 |
| 96 |
|
| 97 |
|
| 98 |
def update(self,dt): |
| 99 |
""" |
| 100 |
sets a new velocity and updates the level set fuction |
| 101 |
|
| 102 |
@param dt: time step forward |
| 103 |
""" |
| 104 |
self.__phi=self.__fc.solve(dt, verbose=False)-self.__diam |
| 105 |
self.__update_count += 1 |
| 106 |
if self.__update_count%self.__reinitialization_each == 0: |
| 107 |
self.__phi=self.__reinitialise(self.__phi) |
| 108 |
self.__fc.setInitialSolution(self.__phi+self.__diam) |
| 109 |
self.__updateInterface() |
| 110 |
|
| 111 |
def setTolerance(self,tolerance=1e-3): |
| 112 |
self.__fc.setTolerance(tolerance) |
| 113 |
|
| 114 |
def __reinitialise(self,phi): |
| 115 |
print "reintialization started:" |
| 116 |
s=self.__makeInterface(phi,1.) |
| 117 |
g=grad(phi) |
| 118 |
w = s*g/(length(g)+EPSILON) |
| 119 |
dtau = 0.5*inf(Function(self.__domain).getSize()) |
| 120 |
print "step size: dt = ",dtau |
| 121 |
iter =0 |
| 122 |
# self.__reinitPde.setValue(q=whereNegative(abs(phi)-self.__h), r=phi) |
| 123 |
# self.__reinitPde.setValue(r=phi) |
| 124 |
while (iter<=self.__reinitialization_steps_max): |
| 125 |
phi_old=phi |
| 126 |
self.__reinitPde.setValue(Y = phi+(dtau/2.)*(s-inner(w,grad(phi)))) |
| 127 |
phi_half = self.__reinitPde.getSolution() |
| 128 |
self.__reinitPde.setValue(Y = phi+dtau*(s-inner(w,grad(phi_half)))) |
| 129 |
phi = self.__reinitPde.getSolution() |
| 130 |
change = Lsup(phi-phi_old)/self.__diam |
| 131 |
print "phi range:",inf(phi), sup(phi) |
| 132 |
print "iteration :", iter, " error:", change |
| 133 |
iter +=1 |
| 134 |
print "reintialization completed." |
| 135 |
return phi |
| 136 |
|
| 137 |
def createParameter(self,value_negative=-1.,value_positive=1): |
| 138 |
out = (value_negative+value_positive)/2. + self.__smoothed_char * ((value_positive-value_negative)/2.) |
| 139 |
return out |
| 140 |
|
| 141 |
|
| 142 |
#================ things from here onwards are not used nor tested: ========================================== |
| 143 |
|
| 144 |
|
| 145 |
def getCharacteristicFunction(self): |
| 146 |
return self.__smoothed_char |
| 147 |
|
| 148 |
|
| 149 |
|
| 150 |
|
| 151 |
def RK2(self,L): |
| 152 |
k0=L(phi) |
| 153 |
phi_1=phi+dt*k0 |
| 154 |
k1=L(phi_1) |
| 155 |
phi=phi+dt/2*(k0+k1) |
| 156 |
def RK3(self,L): |
| 157 |
k0=L(phi) |
| 158 |
phi_1=phi+dt*L(phi) |
| 159 |
k1=L(phi_1) |
| 160 |
phi_2=phi+dt/4*(k0+k1) |
| 161 |
k2=L(phi_2) |
| 162 |
phi=phi+dt/6*(k0+4*k2+K1) |
| 163 |
def TG(self,L): |
| 164 |
k0=L(phi) |
| 165 |
phi_1=phi+dt/2*k0 |
| 166 |
k1=L(phi_1) |
| 167 |
phi=phi+dt*k1 |
| 168 |
|
| 169 |
|
| 170 |
def __reinitialise_old(self): |
| 171 |
#============================================= |
| 172 |
f=0.1 |
| 173 |
f2=0.3 |
| 174 |
h=self.__h |
| 175 |
fs=h.getFunctionSpace() |
| 176 |
vol=self.getVolumeOfNegativeDomain() |
| 177 |
self.__reinitPde.setValue(D=1.0) |
| 178 |
#============ |
| 179 |
grad_phi=grad(self.__phi,fs) |
| 180 |
len_grad_phi=length(grad_phi) |
| 181 |
w=grad_phi/len_grad_phi |
| 182 |
|
| 183 |
# s=wherePositive(self. __makeInterface(2.))-whereNegative(self. __makeInterface(2.)) |
| 184 |
s=self. __makeInterface(2.) |
| 185 |
# s2=whereNegative(len_grad_phi-f2)*(1-(1-len_grad_phi/f2)**2)+(1-whereNegative(len_grad_phi-f2)) |
| 186 |
s2=1. |
| 187 |
s=s*s2 |
| 188 |
dtau=f*inf(h/abs(s)) |
| 189 |
|
| 190 |
#======================== |
| 191 |
diff=1. |
| 192 |
d=1. |
| 193 |
c =0 |
| 194 |
#============================================== |
| 195 |
TVD=integrate(length(grad_phi)) |
| 196 |
print "initial range ",inf(self.__phi),sup(self.__phi),"error:",Lsup(1.-len_grad_phi),"volume =",vol,TVD |
| 197 |
# saveVTK("test.%s.xml"%c,l=length(grad(self.__phi,fs))-1,s=s,phi=self.__phi) |
| 198 |
|
| 199 |
dtau=f*inf(h/abs(s)) |
| 200 |
while c < self.__reinitialization_steps_max: # and abs(diff) >= 0.01: |
| 201 |
# |
| 202 |
grad_phi=grad(self.__phi,fs) |
| 203 |
len_grad_phi=length(grad_phi) |
| 204 |
# |
| 205 |
# s=self.__makeInterface(1.) |
| 206 |
# s2=whereNegative(len_grad_phi-f2)*(1-(1-len_grad_phi/f2)**2)+(1-whereNegative(len_grad_phi-f2)) |
| 207 |
# s=s*s2 |
| 208 |
# phi_on_h=interpolate(self.__phi,fs) |
| 209 |
# self.__reinitPde.setValue(Y =self.__phi+dtau/2*s*(1.-len_grad_phi)) |
| 210 |
# phi_half=self.__reinitPde.getSolution() |
| 211 |
# self.__reinitPde.setValue(Y =self.__phi+dtau*s*(1.-length(grad(phi_half,fs)))) |
| 212 |
# self.__reinitPde.setValue(Y =self.__phi, Y_reduced=dtau*s*(1.-inner(w,grad(phi_half,ReducedFunction(self.__domain))))) |
| 213 |
# self.__reinitPde.setValue(Y =self.__phi+dtau*(s-inner(w,grad_phi)),X=dtau/2*h*(s-inner(w,grad_phi))*w) |
| 214 |
# self.__reinitPde.setValue(Y =self.__phi+dtau*s*(1.-len_grad_phi),X=-dtau/2*h*s**2*grad_phi) |
| 215 |
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) |
| 216 |
|
| 217 |
# 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) |
| 218 |
self.__phi, previous = self.__reinitPde.getSolution(), self.__phi |
| 219 |
# self.__updateInterface() |
| 220 |
|
| 221 |
vol,vol_old=self.getVolumeOfNegativeDomain(),vol |
| 222 |
diff=(vol-vol_old) |
| 223 |
r=Lsup(length(grad(self.__phi))-1.) |
| 224 |
TVD=integrate(length(grad(self.__phi,fs))) |
| 225 |
print "iteration :", c, "range ",inf(self.__phi),sup(self.__phi),"error :",r,"volume change:",diff,TVD |
| 226 |
# saveVTK("test.%s.xml"%(c+1),l=length(grad(self.__phi,fs)),s=s,phi=self.__phi,v=grad(self.__phi,fs)) |
| 227 |
c += 1 |
| 228 |
return |
| 229 |
#============================================== |
| 230 |
f2=0.7 |
| 231 |
h=self.__h |
| 232 |
fs=h.getFunctionSpace() |
| 233 |
vol=self.getVolumeOfNegativeDomain() |
| 234 |
s=abs(self. __makeInterface(2.)) |
| 235 |
grad_phi=grad(self.__phi,fs) |
| 236 |
len_grad_phi=length(grad_phi) |
| 237 |
|
| 238 |
#---- |
| 239 |
# aphi_on_h=abs(interpolate(self.__phi,self.__h.getFunctionSpace())) |
| 240 |
# q=2*self.__h |
| 241 |
# mask_neg = whereNegative(aphi_on_h-q) |
| 242 |
# mask_pos = wherePositive(aphi_on_h-q-self.__h) |
| 243 |
# mask_interface = 1.-mask_neg-mask_pos |
| 244 |
# s2=mask_neg + mask_interface * (q+self.__h-aphi_on_h)/self.__h |
| 245 |
#---- |
| 246 |
m=whereNonPositive(len_grad_phi-f2) |
| 247 |
s2=m*len_grad_phi/f2+(1.-m) |
| 248 |
#---- |
| 249 |
# s2=1. |
| 250 |
#---- |
| 251 |
self.__reinitPde.setValue(D=(1-s)/h**2,Y=(1-s)/h**2*self.__phi) |
| 252 |
self.__reinitPde.setValue(A=(s2*len_grad_phi+(1.-s2))*kronecker(3),X=grad_phi) |
| 253 |
# self.__reinitPde.setValue(A=(len_grad_phi-1)/len_grad_phi*kronecker(3)) |
| 254 |
# self.__reinitPde.setValue(A=kronecker(3), X=grad_phi/len_grad_phi) |
| 255 |
self.__phi = self.__reinitPde.getSolution() |
| 256 |
r=Lsup(length(grad(self.__phi))-1.) |
| 257 |
vol,vol_old=self.getVolumeOfNegativeDomain(),vol |
| 258 |
diff=(vol-vol_old)/vol |
| 259 |
print "iteration :", inf(self.__phi),sup(self.__phi),r,diff |
| 260 |
# saveVTK("test.%s.xml"%0,l=length(grad(self.__phi,fs)),s=s,phi=self.__phi,v=grad(self.__phi,fs),s2=s2) |
| 261 |
return |
| 262 |
#============================================= |
| 263 |
|
| 264 |
#============ |
| 265 |
|
| 266 |
|
| 267 |
def getVolumeOfNegativeDomain(self): |
| 268 |
""" |
| 269 |
return the current volume of domain with phi<0. |
| 270 |
""" |
| 271 |
return integrate((1.-self.__makeInterface(1.))/2.) |
| 272 |
|
| 273 |
|
| 274 |
def getBoundingBoxOfNegativeDomain(self): |
| 275 |
""" |
| 276 |
get the height of the region with phi<0 |
| 277 |
""" |
| 278 |
fs=self.__h.getFunctionSpace() |
| 279 |
mask_phi1=wherePositive(interpolate(self.__phi,fs)) |
| 280 |
mask_phi2=wherePositive(self.__phi) |
| 281 |
x1=fs.getX() |
| 282 |
x2=self.__domain.getX() |
| 283 |
out=[] |
| 284 |
for i in range(fs.getDim()): |
| 285 |
x1_i=x1[i] |
| 286 |
x2_i=x2[i] |
| 287 |
d=2*(sup(x2_i)-inf(x2_i)) |
| 288 |
offset1=d*mask_phi1 |
| 289 |
offset2=d*mask_phi2 |
| 290 |
out.append((min(inf(x1_i+offset1),inf(x2_i+offset2)),max(sup(x1_i-offset1),sup(x2_i-offset2)))) |
| 291 |
return tuple(out) |
| 292 |
|
Parent Directory
| 
Revision Log