/[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 525 - (hide annotations)
Tue Feb 14 06:56:13 2006 UTC (13 years, 6 months ago) by gross
File MIME type: text/x-python
File size: 10437 byte(s)
rob's strange NoPDE class added
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 jgs 149 import escript
13     import linearPDEs
14 jgs 121 import numarray
15 jgs 149 import util
16 jgs 121
17 gross 351 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 gross 396 def getValue(self):
53 gross 409 out=self.__v_mem[0]
54     if len(out)==1:
55     return out[0]
56     else:
57     return out
58    
59 gross 351 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 gross 396
94 gross 351
95 jgs 121 class Projector:
96 jgs 149 """
97     The Projector is a factory which projects a discontiuous function onto a
98     continuous function on the a given domain.
99     """
100 jgs 121 def __init__(self, domain, reduce = True, fast=True):
101     """
102 jgs 149 Create a continuous function space projector for a domain.
103 jgs 121
104 jgs 149 @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 jgs 121 """
108 jgs 149 self.__pde = linearPDEs.LinearPDE(domain)
109 jgs 148 if fast:
110 jgs 149 self.__pde.setSolverMethod(linearPDEs.LinearPDE.LUMPING)
111 jgs 121 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 jgs 149 Projects input_data onto a continuous function
122 jgs 121
123 jgs 149 @param input_data: The input_data to be projected.
124 jgs 121 """
125 gross 525 out=escript.Data(0.,input_data.getShape(),self.__pde.getFunctionSpaceForSolution())
126 jgs 121 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 gross 525 class NoPDE:
154     """
155     solves the following problem for u:
156 jgs 121
157 gross 525 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 jgs 147 class Locator:
267     """
268 jgs 149 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 jgs 147 """
275    
276     def __init__(self,where,x=numarray.zeros((3,))):
277 jgs 149 """
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 jgs 147 self.__function_space=where
284 jgs 121 else:
285 jgs 149 self.__function_space=escript.ContinuousFunction(where)
286 gross 525 self.__id=util.length(self.__function_space.getX()-x[:self.__function_space.getDim()]).mindp()
287 jgs 121
288 jgs 147 def __str__(self):
289 jgs 149 """
290     Returns the coordinates of the Locator as a string.
291     """
292 jgs 147 return "<Locator %s>"%str(self.getX())
293 jgs 121
294 jgs 147 def getFunctionSpace(self):
295 jgs 149 """
296     Returns the function space of the Locator.
297     """
298 jgs 147 return self.__function_space
299    
300 jgs 121 def getId(self):
301 jgs 149 """
302     Returns the identifier of the location.
303     """
304 jgs 121 return self.__id
305    
306 jgs 147 def getX(self):
307 jgs 149 """
308     Returns the exact coordinates of the Locator.
309     """
310 jgs 147 return self(self.getFunctionSpace().getX())
311 jgs 121
312 jgs 147 def __call__(self,data):
313 jgs 149 """
314     Returns the value of data at the Locator of a Data object otherwise
315     the object is returned.
316     """
317 jgs 147 return self.getValue(data)
318 jgs 121
319 jgs 147 def getValue(self,data):
320 jgs 149 """
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 jgs 147 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 jgs 149
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