/[escript]/trunk/escript/py_src/pdetools.py
ViewVC logotype

Contents of /trunk/escript/py_src/pdetools.py

Parent Directory Parent Directory | Revision Log Revision Log


Revision 609 - (show annotations)
Tue Mar 21 09:46:39 2006 UTC (13 years, 7 months ago) by elspeth
File MIME type: text/x-python
File size: 10726 byte(s)
Updated copyright and licence notices.

1 # $Id$
2
3 """
4 Provides some tools related to PDEs.
5
6 Currently includes:
7 - Projector - to project a discontinuous
8 - Locator - to trace values in data objects at a certain location
9 - TimeIntegrationManager - to handel extraplotion in time
10 """
11
12 __copyright__=""" Copyright (c) 2006 by ACcESS MNRF
13 http://www.access.edu.au
14 Primary Business: Queensland, Australia"""
15 __licence__="""Licensed under the Open Software License version 3.0
16 http://www.opensource.org/licences/osl-3.0.php"""
17
18 import escript
19 import linearPDEs
20 import numarray
21 import util
22
23 class TimeIntegrationManager:
24 """
25 a simple mechanism to manage time dependend values.
26
27 typical usage is:
28
29 dt=0.1 # time increment
30 tm=TimeIntegrationManager(inital_value,p=1)
31 while t<1.
32 v_guess=tm.extrapolate(dt) # extrapolate to t+dt
33 v=...
34 tm.checkin(dt,v)
35 t+=dt
36
37 @remark: currently only p=1 is supported.
38 """
39 def __init__(self,*inital_values,**kwargs):
40 """
41 sets up the value manager where inital_value is the initial value and p is order used for extrapolation
42 """
43 if kwargs.has_key("p"):
44 self.__p=kwargs["p"]
45 else:
46 self.__p=1
47 if kwargs.has_key("time"):
48 self.__t=kwargs["time"]
49 else:
50 self.__t=0.
51 self.__v_mem=[inital_values]
52 self.__order=0
53 self.__dt_mem=[]
54 self.__num_val=len(inital_values)
55
56 def getTime(self):
57 return self.__t
58 def getValue(self):
59 out=self.__v_mem[0]
60 if len(out)==1:
61 return out[0]
62 else:
63 return out
64
65 def checkin(self,dt,*values):
66 """
67 adds new values to the manager. the p+1 last value get lost
68 """
69 o=min(self.__order+1,self.__p)
70 self.__order=min(self.__order+1,self.__p)
71 v_mem_new=[values]
72 dt_mem_new=[dt]
73 for i in range(o-1):
74 v_mem_new.append(self.__v_mem[i])
75 dt_mem_new.append(self.__dt_mem[i])
76 v_mem_new.append(self.__v_mem[o-1])
77 self.__order=o
78 self.__v_mem=v_mem_new
79 self.__dt_mem=dt_mem_new
80 self.__t+=dt
81
82 def extrapolate(self,dt):
83 """
84 extrapolates to dt forward in time.
85 """
86 if self.__order==0:
87 out=self.__v_mem[0]
88 else:
89 out=[]
90 for i in range(self.__num_val):
91 out.append((1.+dt/self.__dt_mem[0])*self.__v_mem[0][i]-dt/self.__dt_mem[0]*self.__v_mem[1][i])
92
93 if len(out)==0:
94 return None
95 elif len(out)==1:
96 return out[0]
97 else:
98 return out
99
100
101 class Projector:
102 """
103 The Projector is a factory which projects a discontiuous function onto a
104 continuous function on the a given domain.
105 """
106 def __init__(self, domain, reduce = True, fast=True):
107 """
108 Create a continuous function space projector for a domain.
109
110 @param domain: Domain of the projection.
111 @param reduce: Flag to reduce projection order (default is True)
112 @param fast: Flag to use a fast method based on matrix lumping (default is true)
113 """
114 self.__pde = linearPDEs.LinearPDE(domain)
115 if fast:
116 self.__pde.setSolverMethod(linearPDEs.LinearPDE.LUMPING)
117 self.__pde.setSymmetryOn()
118 self.__pde.setReducedOrderTo(reduce)
119 self.__pde.setValue(D = 1.)
120 return
121
122 def __del__(self):
123 return
124
125 def __call__(self, input_data):
126 """
127 Projects input_data onto a continuous function
128
129 @param input_data: The input_data to be projected.
130 """
131 out=escript.Data(0.,input_data.getShape(),self.__pde.getFunctionSpaceForSolution())
132 if input_data.getRank()==0:
133 self.__pde.setValue(Y = input_data)
134 out=self.__pde.getSolution()
135 elif input_data.getRank()==1:
136 for i0 in range(input_data.getShape()[0]):
137 self.__pde.setValue(Y = input_data[i0])
138 out[i0]=self.__pde.getSolution()
139 elif input_data.getRank()==2:
140 for i0 in range(input_data.getShape()[0]):
141 for i1 in range(input_data.getShape()[1]):
142 self.__pde.setValue(Y = input_data[i0,i1])
143 out[i0,i1]=self.__pde.getSolution()
144 elif input_data.getRank()==3:
145 for i0 in range(input_data.getShape()[0]):
146 for i1 in range(input_data.getShape()[1]):
147 for i2 in range(input_data.getShape()[2]):
148 self.__pde.setValue(Y = input_data[i0,i1,i2])
149 out[i0,i1,i2]=self.__pde.getSolution()
150 else:
151 for i0 in range(input_data.getShape()[0]):
152 for i1 in range(input_data.getShape()[1]):
153 for i2 in range(input_data.getShape()[2]):
154 for i3 in range(input_data.getShape()[3]):
155 self.__pde.setValue(Y = input_data[i0,i1,i2,i3])
156 out[i0,i1,i2,i3]=self.__pde.getSolution()
157 return out
158
159 class NoPDE:
160 """
161 solves the following problem for u:
162
163 M{kronecker[i,j]*D[j]*u[j]=Y[i]}
164
165 with constraint
166
167 M{u[j]=r[j]} where M{q[j]>0}
168
169 where D, Y, r and q are given functions of rank 1.
170
171 In the case of scalars this takes the form
172
173 M{D*u=Y}
174
175 with constraint
176
177 M{u=r} where M{q>0}
178
179 where D, Y, r and q are given scalar functions.
180
181 The constraint is overwriting any other condition.
182
183 @remark: This class is similar to the L{LinearPDE} class with A=B=C=X=0 but has the intention
184 that all input parameter are given in L{Solution} or L{ReducedSolution}. The whole
185 thing is a bit strange and I blame Robert.Woodcock@csiro.au for this.
186 """
187 def __init__(self,domain,D=None,Y=None,q=None,r=None):
188 """
189 initialize the problem
190
191 @param domain: domain of the PDE.
192 @type domain: L{Domain}
193 @param D: coefficient of the solution.
194 @type D: C{float}, C{int}, L{NumArray}, L{Data}
195 @param Y: right hand side
196 @type Y: C{float}, C{int}, L{NumArray}, L{Data}
197 @param q: location of constraints
198 @type q: C{float}, C{int}, L{NumArray}, L{Data}
199 @param r: value of solution at locations of constraints
200 @type r: C{float}, C{int}, L{NumArray}, L{Data}
201 """
202 self.__domain=domain
203 self.__D=D
204 self.__Y=Y
205 self.__q=q
206 self.__r=r
207 self.__u=None
208 self.__function_space=escript.Solution(self.__domain)
209 def setReducedOn(self):
210 """
211 sets the L{FunctionSpace} of the solution to L{ReducedSolution}
212 """
213 self.__function_space=escript.ReducedSolution(self.__domain)
214 self.__u=None
215
216 def setReducedOff(self):
217 """
218 sets the L{FunctionSpace} of the solution to L{Solution}
219 """
220 self.__function_space=escript.Solution(self.__domain)
221 self.__u=None
222
223 def setValue(self,D=None,Y=None,q=None,r=None):
224 """
225 assigns values to the parameters.
226
227 @param D: coefficient of the solution.
228 @type D: C{float}, C{int}, L{NumArray}, L{Data}
229 @param Y: right hand side
230 @type Y: C{float}, C{int}, L{NumArray}, L{Data}
231 @param q: location of constraints
232 @type q: C{float}, C{int}, L{NumArray}, L{Data}
233 @param r: value of solution at locations of constraints
234 @type r: C{float}, C{int}, L{NumArray}, L{Data}
235 """
236 if not D==None:
237 self.__D=D
238 self.__u=None
239 if not Y==None:
240 self.__Y=Y
241 self.__u=None
242 if not q==None:
243 self.__q=q
244 self.__u=None
245 if not r==None:
246 self.__r=r
247 self.__u=None
248
249 def getSolution(self):
250 """
251 returns the solution
252
253 @return: the solution of the problem
254 @rtype: L{Data} object in the L{FunctionSpace} L{Solution} or L{ReducedSolution}.
255 """
256 if self.__u==None:
257 if self.__D==None:
258 raise ValueError,"coefficient D is undefined"
259 D=escript.Data(self.__D,self.__function_space)
260 if D.getRank()>1:
261 raise ValueError,"coefficient D must have rank 0 or 1"
262 if self.__Y==None:
263 self.__u=escript.Data(0.,D.getShape(),self.__function_space)
264 else:
265 self.__u=util.quotient(self.__Y,D)
266 if not self.__q==None:
267 q=util.wherePositive(escript.Data(self.__q,self.__function_space))
268 self.__u*=(1.-q)
269 if not self.__r==None: self.__u+=q*self.__r
270 return self.__u
271
272 class Locator:
273 """
274 Locator provides access to the values of data objects at a given
275 spatial coordinate x.
276
277 In fact, a Locator object finds the sample in the set of samples of a
278 given function space or domain where which is closest to the given
279 point x.
280 """
281
282 def __init__(self,where,x=numarray.zeros((3,))):
283 """
284 Initializes a Locator to access values in Data objects on the Doamin
285 or FunctionSpace where for the sample point which
286 closest to the given point x.
287 """
288 if isinstance(where,escript.FunctionSpace):
289 self.__function_space=where
290 else:
291 self.__function_space=escript.ContinuousFunction(where)
292 self.__id=util.length(self.__function_space.getX()-x[:self.__function_space.getDim()]).mindp()
293
294 def __str__(self):
295 """
296 Returns the coordinates of the Locator as a string.
297 """
298 return "<Locator %s>"%str(self.getX())
299
300 def getFunctionSpace(self):
301 """
302 Returns the function space of the Locator.
303 """
304 return self.__function_space
305
306 def getId(self):
307 """
308 Returns the identifier of the location.
309 """
310 return self.__id
311
312 def getX(self):
313 """
314 Returns the exact coordinates of the Locator.
315 """
316 return self(self.getFunctionSpace().getX())
317
318 def __call__(self,data):
319 """
320 Returns the value of data at the Locator of a Data object otherwise
321 the object is returned.
322 """
323 return self.getValue(data)
324
325 def getValue(self,data):
326 """
327 Returns the value of data at the Locator if data is a Data object
328 otherwise the object is returned.
329 """
330 if isinstance(data,escript.Data):
331 if data.getFunctionSpace()==self.getFunctionSpace():
332 out=data.convertToNumArrayFromDPNo(self.getId()[0],self.getId()[1])
333 else:
334 out=data.interpolate(self.getFunctionSpace()).convertToNumArrayFromDPNo(self.getId()[0],self.getId()[1])
335 if data.getRank()==0:
336 return out[0]
337 else:
338 return out
339 else:
340 return data
341
342 # vim: expandtab shiftwidth=4:

Properties

Name Value
svn:eol-style native
svn:keywords Author Date Id Revision

  ViewVC Help
Powered by ViewVC 1.1.26