/[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 351 - (show annotations)
Tue Dec 13 09:12:15 2005 UTC (13 years, 10 months ago) by gross
File MIME type: text/x-python
File size: 6566 byte(s)
time integration manger 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
53 def checkin(self,dt,*values):
54 """
55 adds new values to the manager. the p+1 last value get lost
56 """
57 o=min(self.__order+1,self.__p)
58 self.__order=min(self.__order+1,self.__p)
59 v_mem_new=[values]
60 dt_mem_new=[dt]
61 for i in range(o-1):
62 v_mem_new.append(self.__v_mem[i])
63 dt_mem_new.append(self.__dt_mem[i])
64 v_mem_new.append(self.__v_mem[o-1])
65 self.__order=o
66 self.__v_mem=v_mem_new
67 self.__dt_mem=dt_mem_new
68 self.__t+=dt
69
70 def extrapolate(self,dt):
71 """
72 extrapolates to dt forward in time.
73 """
74 if self.__order==0:
75 out=self.__v_mem[0]
76 else:
77 out=[]
78 for i in range(self.__num_val):
79 out.append((1.+dt/self.__dt_mem[0])*self.__v_mem[0][i]-dt/self.__dt_mem[0]*self.__v_mem[1][i])
80
81 if len(out)==0:
82 return None
83 elif len(out)==1:
84 return out[0]
85 else:
86 return out
87
88 class Projector:
89 """
90 The Projector is a factory which projects a discontiuous function onto a
91 continuous function on the a given domain.
92 """
93 def __init__(self, domain, reduce = True, fast=True):
94 """
95 Create a continuous function space projector for a domain.
96
97 @param domain: Domain of the projection.
98 @param reduce: Flag to reduce projection order (default is True)
99 @param fast: Flag to use a fast method based on matrix lumping (default is true)
100 """
101 self.__pde = linearPDEs.LinearPDE(domain)
102 if fast:
103 self.__pde.setSolverMethod(linearPDEs.LinearPDE.LUMPING)
104 self.__pde.setSymmetryOn()
105 self.__pde.setReducedOrderTo(reduce)
106 self.__pde.setValue(D = 1.)
107 return
108
109 def __del__(self):
110 return
111
112 def __call__(self, input_data):
113 """
114 Projects input_data onto a continuous function
115
116 @param input_data: The input_data to be projected.
117 """
118 out=escript.Data(0.,input_data.getShape(),what=escript.ContinuousFunction(self.__pde.getDomain()))
119 if input_data.getRank()==0:
120 self.__pde.setValue(Y = input_data)
121 out=self.__pde.getSolution()
122 elif input_data.getRank()==1:
123 for i0 in range(input_data.getShape()[0]):
124 self.__pde.setValue(Y = input_data[i0])
125 out[i0]=self.__pde.getSolution()
126 elif input_data.getRank()==2:
127 for i0 in range(input_data.getShape()[0]):
128 for i1 in range(input_data.getShape()[1]):
129 self.__pde.setValue(Y = input_data[i0,i1])
130 out[i0,i1]=self.__pde.getSolution()
131 elif input_data.getRank()==3:
132 for i0 in range(input_data.getShape()[0]):
133 for i1 in range(input_data.getShape()[1]):
134 for i2 in range(input_data.getShape()[2]):
135 self.__pde.setValue(Y = input_data[i0,i1,i2])
136 out[i0,i1,i2]=self.__pde.getSolution()
137 else:
138 for i0 in range(input_data.getShape()[0]):
139 for i1 in range(input_data.getShape()[1]):
140 for i2 in range(input_data.getShape()[2]):
141 for i3 in range(input_data.getShape()[3]):
142 self.__pde.setValue(Y = input_data[i0,i1,i2,i3])
143 out[i0,i1,i2,i3]=self.__pde.getSolution()
144 return out
145
146
147 class Locator:
148 """
149 Locator provides access to the values of data objects at a given
150 spatial coordinate x.
151
152 In fact, a Locator object finds the sample in the set of samples of a
153 given function space or domain where which is closest to the given
154 point x.
155 """
156
157 def __init__(self,where,x=numarray.zeros((3,))):
158 """
159 Initializes a Locator to access values in Data objects on the Doamin
160 or FunctionSpace where for the sample point which
161 closest to the given point x.
162 """
163 if isinstance(where,escript.FunctionSpace):
164 self.__function_space=where
165 else:
166 self.__function_space=escript.ContinuousFunction(where)
167 self.__id=util.length(x[:self.__function_space.getDim()]-self.__function_space.getX()).mindp()
168
169 def __str__(self):
170 """
171 Returns the coordinates of the Locator as a string.
172 """
173 return "<Locator %s>"%str(self.getX())
174
175 def getFunctionSpace(self):
176 """
177 Returns the function space of the Locator.
178 """
179 return self.__function_space
180
181 def getId(self):
182 """
183 Returns the identifier of the location.
184 """
185 return self.__id
186
187 def getX(self):
188 """
189 Returns the exact coordinates of the Locator.
190 """
191 return self(self.getFunctionSpace().getX())
192
193 def __call__(self,data):
194 """
195 Returns the value of data at the Locator of a Data object otherwise
196 the object is returned.
197 """
198 return self.getValue(data)
199
200 def getValue(self,data):
201 """
202 Returns the value of data at the Locator if data is a Data object
203 otherwise the object is returned.
204 """
205 if isinstance(data,escript.Data):
206 if data.getFunctionSpace()==self.getFunctionSpace():
207 out=data.convertToNumArrayFromDPNo(self.getId()[0],self.getId()[1])
208 else:
209 out=data.interpolate(self.getFunctionSpace()).convertToNumArrayFromDPNo(self.getId()[0],self.getId()[1])
210 if data.getRank()==0:
211 return out[0]
212 else:
213 return out
214 else:
215 return data
216
217 # 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