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

Contents of /temp/escript/py_src/pdetools.py

Parent Directory Parent Directory | Revision Log Revision Log


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