/[escript]/trunk/escriptcore/py_src/levelset.py
ViewVC logotype

Contents of /trunk/escriptcore/py_src/levelset.py

Parent Directory Parent Directory | Revision Log Revision Log


Revision 4821 - (show annotations)
Tue Apr 1 04:58:33 2014 UTC (5 years, 2 months ago) by sshaw
File MIME type: text/x-python
File size: 8299 byte(s)
moved SolverOptions to c++, split into SolverOptions for the options and SolverBuddy as the state as a precursor to per-pde solving... does break some use cases (e.g. pde.getSolverOptions().DIRECT will now fail, new value access is with SolverOptions.DIRECT), examples and documentation updated to match
1
2 ##############################################################################
3 #
4 # Copyright (c) 2003-2014 by University of Queensland
5 # http://www.uq.edu.au
6 #
7 # Primary Business: Queensland, Australia
8 # Licensed under the Open Software License version 3.0
9 # http://www.opensource.org/licenses/osl-3.0.php
10 #
11 # Development until 2012 by Earth Systems Science Computational Center (ESSCC)
12 # Development 2012-2013 by School of Earth Sciences
13 # Development from 2014 by Centre for Geoscience Computing (GeoComp)
14 #
15 ##############################################################################
16
17 __copyright__="""Copyright (c) 2003-2014 by University of Queensland
18 http://www.uq.edu.au
19 Primary Business: Queensland, Australia"""
20 __license__="""Licensed under the Open Software License version 3.0
21 http://www.opensource.org/licenses/osl-3.0.php"""
22 __url__="https://launchpad.net/escript-finley"
23
24 import esys.escriptcore.linearPDEs as lpe
25 import esys.escriptcore.util as es
26 import math
27
28 class LevelSet(object):
29 """
30 The level set method tracking an interface defined by the zero contour of the
31 level set function phi which defines the signed distance of a point x from the
32 interface. The contour phi(x)=0 defines the interface.
33
34 It is assumed that phi(x)<0 defines the volume of interest,
35 phi(x)>0 the outside world.
36 """
37 def __init__(self,phi,reinit_max=10,reinitialize_after=1,smooth=2., useReducedOrder=False):
38 """
39 Sets up the level set method.
40
41 :param phi: the initial level set function
42 :param reinit_max: maximum number of reinitialization steps
43 :param reinitialize_after: ``phi`` is reinitialized every ``reinit_after`` step
44 :param smooth: smoothing width
45 """
46 self.__domain = phi.getDomain()
47 self.__phi = phi
48
49
50 self.__transport=lpe.SingleTransportPDE(self.__domain)
51 if useReducedOrder: self.__transport.setReducedOrderOn()
52 self.__transport.setValue(M=1.0)
53 self.__transport.setInitialSolution(phi)
54
55
56 self.__reinitPDE = lpe.LinearPDE(self.__domain, numEquations=1)
57 self.__reinitPDE.getSolverOptions().setSolverMethod(lpe.SolverOptions.LUMPING)
58 if useReducedOrder: self.__reinitPDE.setReducedOrderOn()
59 self.__reinitPDE.setValue(D=1.0)
60
61 # revise:
62 self.__reinit_max = reinit_max
63 self.__reinit_after = reinitialize_after
64 self.__h = es.inf(self.__domain.getSize())
65 self.__smooth = smooth
66 self.__n_step=0
67
68 def getH(self):
69 """
70 Returns the mesh size.
71 """
72 return self.__h
73
74 def getDomain(self):
75 """
76 Returns the domain.
77 """
78 return self.__domain
79
80 def getAdvectionSolverOptions(self):
81 """
82 Returns the solver options for the interface advective.
83 """
84 return self.__transport.getSolverOptions()
85
86 def getReinitializationSolverOptions(self):
87 """
88 Returns the options of the solver for the reinitialization
89 """
90 return self.__reinitPDE.getSolverOption()
91
92 def getLevelSetFunction(self):
93 """
94 Returns the level set function.
95 """
96 return self.__phi
97
98 def getTimeStepSize(self,flux):
99 """
100 Returns a new ``dt`` for a given ``flux`` using the Courant condition.
101
102 :param flux: flux field
103 """
104 self.__transport.setValue(C=-flux)
105 dt=self.__transport.getSafeTimeStepSize()
106 return dt
107
108 def update(self,dt):
109 """
110 Updates the level set function.
111
112 :param dt: time step forward
113 """
114 self.__phi = self.__transport.getSolution(dt)
115 self.__n_step+=1
116 if self.__n_step%self.__reinit_after ==0:
117 self.__phi = self.__reinitialise(self.__phi)
118 self.__transport.setInitialSolution(self.__phi)
119 return self.__phi
120
121
122 #==============================================================================================
123 def __reinitialise(self, phi):
124 """
125 Reinitializes the level set.
126
127 It solves the PDE...
128
129 :return: reinitialized level set
130 """
131 fs=es.ReducedFunction(self.__domain)
132 dtau = 0.2*self.__h
133 n =0
134 #g=grad(phi,fs)
135 #error = Lsup((1-length(g))*whereNegative(abs(phi.interpolate(fs))-5.0*self.__h))
136 #print(("LevelSet:reinitialization: iteration :", n, " error:", error))
137 #mask = whereNegative(abs(phi)-1.*self.__h)
138 #self.__reinitPDE.setValue(q=mask, r=phi)
139 s= es.sign(phi.interpolate(fs))
140 while (n<=self.__reinit_max):
141 g_phi=es.grad(phi, fs)
142 self.__reinitPDE.setValue(Y = dtau* s * (1 - es.length(g_phi) ), X = - dtau**2/2 * g_phi)
143 phi = phi + self.__reinitPDE.getSolution()
144
145
146 n +=1
147 #g=grad(phi,fs)
148 #error = Lsup((1-length(g))*whereNegative(abs(phi.interpolate(fs))-5.0*self.__h))
149 #print(("LevelSet:reinitialization: iteration :", n, " error:", error))
150 return phi
151
152
153 def getVolume(self):
154 """
155 Returns the volume of the *phi(x)<0* region.
156 """
157 return integrate(es.whereNegative(self.__phi.interpolate(es.Function(self.__domain))))
158
159
160 def getJumpingParameter(self, param_neg=-1, param_pos=1, phi=None):
161 """
162 Creates a function with ``param_neg`` where ``phi<0`` and ``param_pos``
163 where ``phi>0`` (no smoothing).
164
165 :param param_neg: value of parameter on the negative side (phi<0)
166 :param param_pos: value of parameter on the positive side (phi>0)
167 :param phi: level set function to be used. If not present the current
168 level set is used.
169 """
170 mask_neg = es.whereNegative(self.__phi)
171 mask_pos = es.whereNonNegative(self.__phi)
172 param = param_pos*mask_pos + param_neg*mask_neg
173 return param
174
175 def getSmoothedParameter(self, param_neg=-1, param_pos=1, phi=None, smoothing_width=None):
176 """
177 Creates a smoothed function with ``param_neg`` where ``phi<0`` and
178 ``param_pos`` where ``phi>0`` which is smoothed over a length
179 ``smoothing_width`` across the interface.
180
181 :param smoothing_width: width of the smoothing zone relative to mesh size.
182 If not present the initial value of ``smooth`` is
183 used.
184 """
185 if smoothing_width==None: smoothing_width = self.__smooth
186 if phi==None: phi=self.__phi
187 s=self.getSmoothedJump(phi,smoothing_width)
188 return ((param_pos-param_neg)*s+param_pos+param_neg)/2
189
190 def getSmoothedJump(self,phi=None,smoothing_width=None):
191 """
192 Creates a smooth interface from -1 to 1 over the length
193 *2*h*smoothing_width* where -1 is used where the level set is negative
194 and 1 where the level set is 1.
195 """
196 if smoothing_width==None: smoothing_width = self.__smooth
197 if phi==None: phi = self.__phi
198 s=smoothing_width*self.__h
199 phi_on_h=es.interpolate(phi,es.Function(self.__domain))
200 mask_neg = es.whereNonNegative(-s-phi_on_h)
201 mask_pos = es.whereNonNegative(phi_on_h-s)
202 mask_interface = 1.-mask_neg-mask_pos
203 interface=phi_on_h/s
204 return - mask_neg + mask_pos + mask_interface * interface
205
206 def getInterface(self,phi=None,smoothing_width=None):
207 """
208 creates a characteristic function which is 1 over the over the length
209 *2*h*smoothing_width* around the interface and zero elsewhere
210 """
211 if smoothing_width==None: smoothing_width = self.__smooth
212 if phi==None: phi = self.__phi
213 s=smoothing_width*self.__h
214 phi_on_h=es.interpolate(phi,es.Function(self.__domain))
215 return es.whereNegative(abs(phi_on_h)-s)
216
217 def makeCharacteristicFunction(self, contour=0, phi=None, positiveSide=True, smoothing_width=None):
218 """
219 Makes a smooth characteristic function of the region ``phi(x)>contour`` if
220 ``positiveSide`` and ``phi(x)<contour`` otherwise.
221
222 :param phi: level set function to be used. If not present the current
223 level set is used.
224 :param smoothing_width: width of the smoothing zone relative to mesh size.
225 If not present the initial value of ``smooth`` is
226 used.
227 """
228 if phi==None: phi=self.__phi
229 if smoothing_width == None: smoothing_width=self.__smooth
230 s=self.getSmoothedJump(phi=phi-contour,smoothing_width=smoothing_width)
231 if positiveSide:
232 return (1+s)/2
233 else:
234 return (1-s)/2

  ViewVC Help
Powered by ViewVC 1.1.26