/[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 782 - (hide annotations)
Tue Jul 18 00:47:47 2006 UTC (13 years, 1 month ago) by bcumming
File MIME type: text/x-python
File size: 11391 byte(s)
Large number of changes to Finley for meshing in MPI.

- optimisation and neatening up of rectcanglular mesh generation code
- first and second order 1D, 2D and 3D rectangular meshes are now
  available in finley and escript using MPI.
- reduced meshes now generated in MPI, and interpolation to and from 
  reduced data types now supported.  

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