/[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 396 - (hide annotations)
Wed Dec 21 05:08:25 2005 UTC (13 years, 11 months ago) by gross
File MIME type: text/x-python
File size: 6618 byte(s)
tests for clip, maximum, minimum added
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     return self.__v_mem[0]
54 gross 351 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 gross 396
89 gross 351
90 jgs 121 class Projector:
91 jgs 149 """
92     The Projector is a factory which projects a discontiuous function onto a
93     continuous function on the a given domain.
94     """
95 jgs 121 def __init__(self, domain, reduce = True, fast=True):
96     """
97 jgs 149 Create a continuous function space projector for a domain.
98 jgs 121
99 jgs 149 @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 jgs 121 """
103 jgs 149 self.__pde = linearPDEs.LinearPDE(domain)
104 jgs 148 if fast:
105 jgs 149 self.__pde.setSolverMethod(linearPDEs.LinearPDE.LUMPING)
106 jgs 121 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 jgs 149 Projects input_data onto a continuous function
117 jgs 121
118 jgs 149 @param input_data: The input_data to be projected.
119 jgs 121 """
120 jgs 149 out=escript.Data(0.,input_data.getShape(),what=escript.ContinuousFunction(self.__pde.getDomain()))
121 jgs 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 jgs 147 class Locator:
150     """
151 jgs 149 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 jgs 147 """
158    
159     def __init__(self,where,x=numarray.zeros((3,))):
160 jgs 149 """
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 jgs 147 self.__function_space=where
167 jgs 121 else:
168 jgs 149 self.__function_space=escript.ContinuousFunction(where)
169     self.__id=util.length(x[:self.__function_space.getDim()]-self.__function_space.getX()).mindp()
170 jgs 121
171 jgs 147 def __str__(self):
172 jgs 149 """
173     Returns the coordinates of the Locator as a string.
174     """
175 jgs 147 return "<Locator %s>"%str(self.getX())
176 jgs 121
177 jgs 147 def getFunctionSpace(self):
178 jgs 149 """
179     Returns the function space of the Locator.
180     """
181 jgs 147 return self.__function_space
182    
183 jgs 121 def getId(self):
184 jgs 149 """
185     Returns the identifier of the location.
186     """
187 jgs 121 return self.__id
188    
189 jgs 147 def getX(self):
190 jgs 149 """
191     Returns the exact coordinates of the Locator.
192     """
193 jgs 147 return self(self.getFunctionSpace().getX())
194 jgs 121
195 jgs 147 def __call__(self,data):
196 jgs 149 """
197     Returns the value of data at the Locator of a Data object otherwise
198     the object is returned.
199     """
200 jgs 147 return self.getValue(data)
201 jgs 121
202 jgs 147 def getValue(self,data):
203 jgs 149 """
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 jgs 147 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 jgs 149
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