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

Diff of /trunk/escript/py_src/pdetools.py

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

trunk/esys2/escript/py_src/pdetools.py revision 148 by jgs, Tue Aug 23 01:24:31 2005 UTC trunk/escript/py_src/pdetools.py revision 782 by bcumming, Tue Jul 18 00:47:47 2006 UTC
# Line 1  Line 1 
1  # $Id$  # $Id$
2    
3  """  """
4  provides a some tools related to PDEs currently includes:  Provides some tools related to PDEs.
5    
6       Projector - to project a discontinuous  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    @var __author__: name of author
12    @var __copyright__: copyrights
13    @var __license__: licence agreement
14    @var __url__: url entry point on documentation
15    @var __version__: version
16    @var __date__: date of the version
17    """
18    
19    __author__="Lutz Gross, l.gross@uq.edu.au"
20    __copyright__="""  Copyright (c) 2006 by ACcESS MNRF
21                        http://www.access.edu.au
22                    Primary Business: Queensland, Australia"""
23    __license__="""Licensed under the Open Software License version 3.0
24                 http://www.opensource.org/licenses/osl-3.0.php"""
25    __url__="http://www.iservo.edu.au/esys"
26    __version__="$Revision$"
27    __date__="$Date$"
28    
 """  
29    
30  from escript import *  import escript
31  from linearPDEs import LinearPDE  import linearPDEs
32  import numarray  import numarray
33    import util
34    
35    class TimeIntegrationManager:
36      """
37      a simple mechanism to manage time dependend values.
38    
39      typical usage is::
40    
41         dt=0.1 # time increment
42         tm=TimeIntegrationManager(inital_value,p=1)
43         while t<1.
44             v_guess=tm.extrapolate(dt) # extrapolate to t+dt
45             v=...
46             tm.checkin(dt,v)
47             t+=dt
48    
49      @note: currently only p=1 is supported.
50      """
51      def __init__(self,*inital_values,**kwargs):
52         """
53         sets up the value manager where inital_value is the initial value and p is order used for extrapolation
54         """
55         if kwargs.has_key("p"):
56                self.__p=kwargs["p"]
57         else:
58                self.__p=1
59         if kwargs.has_key("time"):
60                self.__t=kwargs["time"]
61         else:
62                self.__t=0.
63         self.__v_mem=[inital_values]
64         self.__order=0
65         self.__dt_mem=[]
66         self.__num_val=len(inital_values)
67    
68      def getTime(self):
69          return self.__t
70      def getValue(self):
71          out=self.__v_mem[0]
72          if len(out)==1:
73              return out[0]
74          else:
75              return out
76    
77      def checkin(self,dt,*values):
78          """
79          adds new values to the manager. the p+1 last value get lost
80          """
81          o=min(self.__order+1,self.__p)
82          self.__order=min(self.__order+1,self.__p)
83          v_mem_new=[values]
84          dt_mem_new=[dt]
85          for i in range(o-1):
86             v_mem_new.append(self.__v_mem[i])
87             dt_mem_new.append(self.__dt_mem[i])
88          v_mem_new.append(self.__v_mem[o-1])
89          self.__order=o
90          self.__v_mem=v_mem_new
91          self.__dt_mem=dt_mem_new
92          self.__t+=dt
93    
94      def extrapolate(self,dt):
95          """
96          extrapolates to dt forward in time.
97          """
98          if self.__order==0:
99             out=self.__v_mem[0]
100          else:
101            out=[]
102            for i in range(self.__num_val):
103               out.append((1.+dt/self.__dt_mem[0])*self.__v_mem[0][i]-dt/self.__dt_mem[0]*self.__v_mem[1][i])
104    
105          if len(out)==0:
106             return None
107          elif len(out)==1:
108             return out[0]
109          else:
110             return out
111    
112    
113  class Projector:  class Projector:
114    """The Projector is a factory which projects a discontiuous function onto a continuous function on the a given domain"""    """
115      The Projector is a factory which projects a discontiuous function onto a
116      continuous function on the a given domain.
117      """
118    def __init__(self, domain, reduce = True, fast=True):    def __init__(self, domain, reduce = True, fast=True):
119      """      """
120      @brief Create a continuous function space projector for a domain.      Create a continuous function space projector for a domain.
121    
122      @param domain Domain of the projection.      @param domain: Domain of the projection.
123      @param reduce Flag to reduce projection order (default is True)      @param reduce: Flag to reduce projection order (default is True)
124      @param fast Flag to use a fast method based on matrix lumping (default is true)      @param fast: Flag to use a fast method based on matrix lumping (default is true)
125      """      """
126      self.__pde = LinearPDE(domain)      self.__pde = linearPDEs.LinearPDE(domain)
127      if fast:      if fast:
128        self.__pde.setSolverMethod(LinearPDE.LUMPING)        self.__pde.setSolverMethod(linearPDEs.LinearPDE.LUMPING)
129      self.__pde.setSymmetryOn()      self.__pde.setSymmetryOn()
130      self.__pde.setReducedOrderTo(reduce)      self.__pde.setReducedOrderTo(reduce)
131      self.__pde.setValue(D = 1.)      self.__pde.setValue(D = 1.)
# Line 35  class Projector: Line 136  class Projector:
136    
137    def __call__(self, input_data):    def __call__(self, input_data):
138      """      """
139      @brief projects input_data onto a continuous function      Projects input_data onto a continuous function
140    
141      @param input_data  The input_data to be projected.      @param input_data: The input_data to be projected.
142      """      """
143      out=Data(0.,input_data.getShape(),what=ContinuousFunction(self.__pde.getDomain()))      out=escript.Data(0.,input_data.getShape(),self.__pde.getFunctionSpaceForSolution())
144      if input_data.getRank()==0:      if input_data.getRank()==0:
145          self.__pde.setValue(Y = input_data)          self.__pde.setValue(Y = input_data)
146          out=self.__pde.getSolution()          out=self.__pde.getSolution()
# Line 67  class Projector: Line 168  class Projector:
168                      out[i0,i1,i2,i3]=self.__pde.getSolution()                      out[i0,i1,i2,i3]=self.__pde.getSolution()
169      return out      return out
170    
171    class NoPDE:
 class Locator:  
172       """       """
173         solves the following problem for u:
174    
175         M{kronecker[i,j]*D[j]*u[j]=Y[i]}
176    
177         with constraint
178    
179         M{u[j]=r[j]}  where M{q[j]>0}
180    
181          Locator provides a  access the values of data objects at a given spatial coordinate x.       where D, Y, r and q are given functions of rank 1.
         In fact, a Locator object finds the sample in the set of samples of a given function space or domain where which is closest  
         to the given point x.  
182    
183         In the case of scalars this takes the form
184    
185         M{D*u=Y}
186    
187         with constraint
188    
189         M{u=r}  where M{q>0}
190    
191         where D, Y, r and q are given scalar functions.
192    
193         The constraint is overwriting any other condition.
194    
195         @note: This class is similar to the L{linearPDEs.LinearPDE} class with A=B=C=X=0 but has the intention
196                that all input parameter are given in L{Solution} or L{ReducedSolution}. The whole
197                thing is a bit strange and I blame Robert.Woodcock@csiro.au for this.
198         """
199         def __init__(self,domain,D=None,Y=None,q=None,r=None):
200             """
201             initialize the problem
202    
203             @param domain: domain of the PDE.
204             @type domain: L{Domain}
205             @param D: coefficient of the solution.
206             @type D: C{float}, C{int}, L{numarray.NumArray}, L{Data}
207             @param Y: right hand side
208             @type Y: C{float}, C{int}, L{numarray.NumArray}, L{Data}
209             @param q: location of constraints
210             @type q: C{float}, C{int}, L{numarray.NumArray}, L{Data}
211             @param r: value of solution at locations of constraints
212             @type r: C{float}, C{int}, L{numarray.NumArray}, L{Data}
213             """
214             self.__domain=domain
215             self.__D=D
216             self.__Y=Y
217             self.__q=q
218             self.__r=r
219             self.__u=None
220             self.__function_space=escript.Solution(self.__domain)
221         def setReducedOn(self):
222             """
223             sets the L{FunctionSpace} of the solution to L{ReducedSolution}
224             """
225             self.__function_space=escript.ReducedSolution(self.__domain)
226             self.__u=None
227    
228         def setReducedOff(self):
229             """
230             sets the L{FunctionSpace} of the solution to L{Solution}
231             """
232             self.__function_space=escript.Solution(self.__domain)
233             self.__u=None
234            
235         def setValue(self,D=None,Y=None,q=None,r=None):
236             """
237             assigns values to the parameters.
238    
239             @param D: coefficient of the solution.
240             @type D: C{float}, C{int}, L{numarray.NumArray}, L{Data}
241             @param Y: right hand side
242             @type Y: C{float}, C{int}, L{numarray.NumArray}, L{Data}
243             @param q: location of constraints
244             @type q: C{float}, C{int}, L{numarray.NumArray}, L{Data}
245             @param r: value of solution at locations of constraints
246             @type r: C{float}, C{int}, L{numarray.NumArray}, L{Data}
247             """
248             if not D==None:
249                self.__D=D
250                self.__u=None
251             if not Y==None:
252                self.__Y=Y
253                self.__u=None
254             if not q==None:
255                self.__q=q
256                self.__u=None
257             if not r==None:
258                self.__r=r
259                self.__u=None
260    
261         def getSolution(self):
262             """
263             returns the solution
264            
265             @return: the solution of the problem
266             @rtype: L{Data} object in the L{FunctionSpace} L{Solution} or L{ReducedSolution}.
267             """
268             if self.__u==None:
269                if self.__D==None:
270                   raise ValueError,"coefficient D is undefined"
271                D=escript.Data(self.__D,self.__function_space)
272                if D.getRank()>1:
273                   raise ValueError,"coefficient D must have rank 0 or 1"
274                if self.__Y==None:
275                   self.__u=escript.Data(0.,D.getShape(),self.__function_space)
276                else:
277                   self.__u=util.quotient(self.__Y,D)
278                if not self.__q==None:
279                    q=util.wherePositive(escript.Data(self.__q,self.__function_space))
280                    self.__u*=(1.-q)
281                    if not self.__r==None: self.__u+=q*self.__r
282             return self.__u
283                
284    class Locator:
285         """
286         Locator provides access to the values of data objects at a given
287         spatial coordinate x.  
288        
289         In fact, a Locator object finds the sample in the set of samples of a
290         given function space or domain where which is closest to the given
291         point x.
292       """       """
293    
294       def __init__(self,where,x=numarray.zeros((3,))):       def __init__(self,where,x=numarray.zeros((3,))):
295         """initializes a Locator to access values in Data objects on the Doamin or FunctionSpace where for the sample point which         """
296            closest to the given point x"""         Initializes a Locator to access values in Data objects on the Doamin
297         if isinstance(where,FunctionSpace):         or FunctionSpace where for the sample point which
298           closest to the given point x.
299           """
300           if isinstance(where,escript.FunctionSpace):
301            self.__function_space=where            self.__function_space=where
302         else:         else:
303            self.__function_space=ContinuousFunction(where)            self.__function_space=escript.ContinuousFunction(where)
304         self.__id=length(x[:self.__function_space.getDim()]-self.__function_space.getX()).mindp()         self.__id=util.length(self.__function_space.getX()-x[:self.__function_space.getDim()]).mindp()
305    
306       def __str__(self):       def __str__(self):
307         """returns the coordinates of the Locator as a string"""         """
308           Returns the coordinates of the Locator as a string.
309           """
310         return "<Locator %s>"%str(self.getX())         return "<Locator %s>"%str(self.getX())
311    
312       def getFunctionSpace(self):       def getFunctionSpace(self):
313          """returns the function space of the Locator"""          """
314        Returns the function space of the Locator.
315        """
316          return self.__function_space          return self.__function_space
317    
318       def getId(self):       def getId(self):
319          """returns the identifier of the location"""          """
320        Returns the identifier of the location.
321        """
322          return self.__id          return self.__id
323    
324       def getX(self):       def getX(self):
325          """returns the exact coordinates of the Locator"""          """
326        Returns the exact coordinates of the Locator.
327        """
328          return self(self.getFunctionSpace().getX())          return self(self.getFunctionSpace().getX())
329    
330       def __call__(self,data):       def __call__(self,data):
331          """returns the value of data at the Locator of a Data object otherwise the object is returned."""          """
332        Returns the value of data at the Locator of a Data object otherwise
333        the object is returned.
334        """
335          return self.getValue(data)          return self.getValue(data)
336    
337       def getValue(self,data):       def getValue(self,data):
338          """returns the value of data at the Locator if data is a Data object otherwise the object is returned."""          """
339          if isinstance(data,Data):      Returns the value of data at the Locator if data is a Data object
340        otherwise the object is returned.
341        """
342            if isinstance(data,escript.Data):
343             if data.getFunctionSpace()==self.getFunctionSpace():             if data.getFunctionSpace()==self.getFunctionSpace():
344               out=data.convertToNumArrayFromDPNo(self.getId()[0],self.getId()[1])               out=data.convertToNumArrayFromDPNo(self.getId()[0],self.getId()[1])
345                 #out=data.convertToNumArrayFromDPNo(self.getId()[0],self.getId()[1],self.getId()[2])
346             else:             else:
347               out=data.interpolate(self.getFunctionSpace()).convertToNumArrayFromDPNo(self.getId()[0],self.getId()[1])               out=data.interpolate(self.getFunctionSpace()).convertToNumArrayFromDPNo(self.getId()[0],self.getId()[1])
348                 #out=data.interpolate(self.getFunctionSpace()).convertToNumArrayFromDPNo(self.getId()[0],self.getId()[1],self.getId()[2])
349             if data.getRank()==0:             if data.getRank()==0:
350                return out[0]                return out[0]
351             else:             else:
352                return out                return out
353          else:          else:
354             return data             return data
355    
356    # vim: expandtab shiftwidth=4:

Legend:
Removed from v.148  
changed lines
  Added in v.782

  ViewVC Help
Powered by ViewVC 1.1.26