/[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 409 - (show annotations)
Fri Dec 23 01:13:41 2005 UTC (13 years, 7 months ago) by gross
File MIME type: text/x-python
File size: 6695 byte(s)
bug fixed
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(),what=escript.ContinuousFunction(self.__pde.getDomain()))
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
154 class Locator:
155 """
156 Locator provides access to the values of data objects at a given
157 spatial coordinate x.
158
159 In fact, a Locator object finds the sample in the set of samples of a
160 given function space or domain where which is closest to the given
161 point x.
162 """
163
164 def __init__(self,where,x=numarray.zeros((3,))):
165 """
166 Initializes a Locator to access values in Data objects on the Doamin
167 or FunctionSpace where for the sample point which
168 closest to the given point x.
169 """
170 if isinstance(where,escript.FunctionSpace):
171 self.__function_space=where
172 else:
173 self.__function_space=escript.ContinuousFunction(where)
174 self.__id=util.length(x[:self.__function_space.getDim()]-self.__function_space.getX()).mindp()
175
176 def __str__(self):
177 """
178 Returns the coordinates of the Locator as a string.
179 """
180 return "<Locator %s>"%str(self.getX())
181
182 def getFunctionSpace(self):
183 """
184 Returns the function space of the Locator.
185 """
186 return self.__function_space
187
188 def getId(self):
189 """
190 Returns the identifier of the location.
191 """
192 return self.__id
193
194 def getX(self):
195 """
196 Returns the exact coordinates of the Locator.
197 """
198 return self(self.getFunctionSpace().getX())
199
200 def __call__(self,data):
201 """
202 Returns the value of data at the Locator of a Data object otherwise
203 the object is returned.
204 """
205 return self.getValue(data)
206
207 def getValue(self,data):
208 """
209 Returns the value of data at the Locator if data is a Data object
210 otherwise the object is returned.
211 """
212 if isinstance(data,escript.Data):
213 if data.getFunctionSpace()==self.getFunctionSpace():
214 out=data.convertToNumArrayFromDPNo(self.getId()[0],self.getId()[1])
215 else:
216 out=data.interpolate(self.getFunctionSpace()).convertToNumArrayFromDPNo(self.getId()[0],self.getId()[1])
217 if data.getRank()==0:
218 return out[0]
219 else:
220 return out
221 else:
222 return data
223
224 # 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