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