/[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 409 - (hide annotations)
Fri Dec 23 01:13:41 2005 UTC (13 years, 10 months ago) by gross
File MIME type: text/x-python
File size: 6695 byte(s)
bug fixed
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 jgs 149 out=escript.Data(0.,input_data.getShape(),what=escript.ContinuousFunction(self.__pde.getDomain()))
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    
154 jgs 147 class Locator:
155     """
156 jgs 149 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 jgs 147 """
163    
164     def __init__(self,where,x=numarray.zeros((3,))):
165 jgs 149 """
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 jgs 147 self.__function_space=where
172 jgs 121 else:
173 jgs 149 self.__function_space=escript.ContinuousFunction(where)
174     self.__id=util.length(x[:self.__function_space.getDim()]-self.__function_space.getX()).mindp()
175 jgs 121
176 jgs 147 def __str__(self):
177 jgs 149 """
178     Returns the coordinates of the Locator as a string.
179     """
180 jgs 147 return "<Locator %s>"%str(self.getX())
181 jgs 121
182 jgs 147 def getFunctionSpace(self):
183 jgs 149 """
184     Returns the function space of the Locator.
185     """
186 jgs 147 return self.__function_space
187    
188 jgs 121 def getId(self):
189 jgs 149 """
190     Returns the identifier of the location.
191     """
192 jgs 121 return self.__id
193    
194 jgs 147 def getX(self):
195 jgs 149 """
196     Returns the exact coordinates of the Locator.
197     """
198 jgs 147 return self(self.getFunctionSpace().getX())
199 jgs 121
200 jgs 147 def __call__(self,data):
201 jgs 149 """
202     Returns the value of data at the Locator of a Data object otherwise
203     the object is returned.
204     """
205 jgs 147 return self.getValue(data)
206 jgs 121
207 jgs 147 def getValue(self,data):
208 jgs 149 """
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 jgs 147 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 jgs 149
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