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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 614 - (hide annotations)
Wed Mar 22 01:37:07 2006 UTC (13 years, 4 months ago) by elspeth
File MIME type: text/x-python
File size: 10726 byte(s)
Corrected spelling of 'license' in url so that the link actually points to the license.

1 jgs 121 # $Id$
2    
3     """
4 jgs 149 Provides some tools related to PDEs.
5 jgs 121
6 jgs 149 Currently includes:
7     - Projector - to project a discontinuous
8 gross 351 - Locator - to trace values in data objects at a certain location
9     - TimeIntegrationManager - to handel extraplotion in time
10 jgs 121 """
11    
12 elspeth 609 __copyright__=""" Copyright (c) 2006 by ACcESS MNRF
13     http://www.access.edu.au
14     Primary Business: Queensland, Australia"""
15 elspeth 614 __license__="""Licensed under the Open Software License version 3.0
16     http://www.opensource.org/licenses/osl-3.0.php"""
17 elspeth 609
18 jgs 149 import escript
19     import linearPDEs
20 jgs 121 import numarray
21 jgs 149 import util
22 jgs 121
23 gross 351 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 gross 396 def getValue(self):
59 gross 409 out=self.__v_mem[0]
60     if len(out)==1:
61     return out[0]
62     else:
63     return out
64    
65 gross 351 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 gross 396
100 gross 351
101 jgs 121 class Projector:
102 jgs 149 """
103     The Projector is a factory which projects a discontiuous function onto a
104     continuous function on the a given domain.
105     """
106 jgs 121 def __init__(self, domain, reduce = True, fast=True):
107     """
108 jgs 149 Create a continuous function space projector for a domain.
109 jgs 121
110 jgs 149 @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 jgs 121 """
114 jgs 149 self.__pde = linearPDEs.LinearPDE(domain)
115 jgs 148 if fast:
116 jgs 149 self.__pde.setSolverMethod(linearPDEs.LinearPDE.LUMPING)
117 jgs 121 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 jgs 149 Projects input_data onto a continuous function
128 jgs 121
129 jgs 149 @param input_data: The input_data to be projected.
130 jgs 121 """
131 gross 525 out=escript.Data(0.,input_data.getShape(),self.__pde.getFunctionSpaceForSolution())
132 jgs 121 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 gross 525 class NoPDE:
160     """
161     solves the following problem for u:
162 jgs 121
163 gross 525 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 jgs 147 class Locator:
273     """
274 jgs 149 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 jgs 147 """
281    
282     def __init__(self,where,x=numarray.zeros((3,))):
283 jgs 149 """
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 jgs 147 self.__function_space=where
290 jgs 121 else:
291 jgs 149 self.__function_space=escript.ContinuousFunction(where)
292 gross 525 self.__id=util.length(self.__function_space.getX()-x[:self.__function_space.getDim()]).mindp()
293 jgs 121
294 jgs 147 def __str__(self):
295 jgs 149 """
296     Returns the coordinates of the Locator as a string.
297     """
298 jgs 147 return "<Locator %s>"%str(self.getX())
299 jgs 121
300 jgs 147 def getFunctionSpace(self):
301 jgs 149 """
302     Returns the function space of the Locator.
303     """
304 jgs 147 return self.__function_space
305    
306 jgs 121 def getId(self):
307 jgs 149 """
308     Returns the identifier of the location.
309     """
310 jgs 121 return self.__id
311    
312 jgs 147 def getX(self):
313 jgs 149 """
314     Returns the exact coordinates of the Locator.
315     """
316 jgs 147 return self(self.getFunctionSpace().getX())
317 jgs 121
318 jgs 147 def __call__(self,data):
319 jgs 149 """
320     Returns the value of data at the Locator of a Data object otherwise
321     the object is returned.
322     """
323 jgs 147 return self.getValue(data)
324 jgs 121
325 jgs 147 def getValue(self,data):
326 jgs 149 """
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 jgs 147 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 jgs 149
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