/[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 149 by jgs, Thu Sep 1 03:31:39 2005 UTC trunk/escript/py_src/pdetools.py revision 609 by elspeth, Tue Mar 21 09:46:39 2006 UTC
# Line 5  Provides some tools related to PDEs. Line 5  Provides some tools related to PDEs.
5    
6  Currently includes:  Currently includes:
7      - Projector - to project a discontinuous      - 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    __copyright__="""  Copyright (c) 2006 by ACcESS MNRF
13                        http://www.access.edu.au
14                    Primary Business: Queensland, Australia"""
15    __licence__="""Licensed under the Open Software License version 3.0
16                 http://www.opensource.org/licences/osl-3.0.php"""
17    
18  import escript  import escript
19  import linearPDEs  import linearPDEs
20  import numarray  import numarray
21  import util  import util
22    
23    class TimeIntegrationManager:
24      """
25      a simple mechanism to manage time dependend values.
26    
27      typical usage is:
28    
29      dt=0.1 # time increment
30      tm=TimeIntegrationManager(inital_value,p=1)
31      while t<1.
32          v_guess=tm.extrapolate(dt) # extrapolate to t+dt
33          v=...
34          tm.checkin(dt,v)
35          t+=dt
36    
37      @remark: currently only p=1 is supported.
38      """
39      def __init__(self,*inital_values,**kwargs):
40         """
41         sets up the value manager where inital_value is the initial value and p is order used for extrapolation
42         """
43         if kwargs.has_key("p"):
44                self.__p=kwargs["p"]
45         else:
46                self.__p=1
47         if kwargs.has_key("time"):
48                self.__t=kwargs["time"]
49         else:
50                self.__t=0.
51         self.__v_mem=[inital_values]
52         self.__order=0
53         self.__dt_mem=[]
54         self.__num_val=len(inital_values)
55    
56      def getTime(self):
57          return self.__t
58      def getValue(self):
59          out=self.__v_mem[0]
60          if len(out)==1:
61              return out[0]
62          else:
63              return out
64    
65      def checkin(self,dt,*values):
66          """
67          adds new values to the manager. the p+1 last value get lost
68          """
69          o=min(self.__order+1,self.__p)
70          self.__order=min(self.__order+1,self.__p)
71          v_mem_new=[values]
72          dt_mem_new=[dt]
73          for i in range(o-1):
74             v_mem_new.append(self.__v_mem[i])
75             dt_mem_new.append(self.__dt_mem[i])
76          v_mem_new.append(self.__v_mem[o-1])
77          self.__order=o
78          self.__v_mem=v_mem_new
79          self.__dt_mem=dt_mem_new
80          self.__t+=dt
81    
82      def extrapolate(self,dt):
83          """
84          extrapolates to dt forward in time.
85          """
86          if self.__order==0:
87             out=self.__v_mem[0]
88          else:
89            out=[]
90            for i in range(self.__num_val):
91               out.append((1.+dt/self.__dt_mem[0])*self.__v_mem[0][i]-dt/self.__dt_mem[0]*self.__v_mem[1][i])
92    
93          if len(out)==0:
94             return None
95          elif len(out)==1:
96             return out[0]
97          else:
98             return out
99    
100    
101  class Projector:  class Projector:
102    """    """
103    The Projector is a factory which projects a discontiuous function onto a    The Projector is a factory which projects a discontiuous function onto a
# Line 42  class Projector: Line 128  class Projector:
128    
129      @param input_data: The input_data to be projected.      @param input_data: The input_data to be projected.
130      """      """
131      out=escript.Data(0.,input_data.getShape(),what=escript.ContinuousFunction(self.__pde.getDomain()))      out=escript.Data(0.,input_data.getShape(),self.__pde.getFunctionSpaceForSolution())
132      if input_data.getRank()==0:      if input_data.getRank()==0:
133          self.__pde.setValue(Y = input_data)          self.__pde.setValue(Y = input_data)
134          out=self.__pde.getSolution()          out=self.__pde.getSolution()
# Line 70  class Projector: Line 156  class Projector:
156                      out[i0,i1,i2,i3]=self.__pde.getSolution()                      out[i0,i1,i2,i3]=self.__pde.getSolution()
157      return out      return out
158    
159    class NoPDE:
160         """
161         solves the following problem for u:
162    
163         M{kronecker[i,j]*D[j]*u[j]=Y[i]}
164    
165         with constraint
166    
167         M{u[j]=r[j]}  where M{q[j]>0}
168    
169         where D, Y, r and q are given functions of rank 1.
170    
171         In the case of scalars this takes the form
172    
173         M{D*u=Y}
174    
175         with constraint
176    
177         M{u=r}  where M{q>0}
178    
179         where D, Y, r and q are given scalar functions.
180    
181         The constraint is overwriting any other condition.
182    
183         @remark: This class is similar to the L{LinearPDE} class with A=B=C=X=0 but has the intention
184                  that all input parameter are given in L{Solution} or L{ReducedSolution}. The whole
185                  thing is a bit strange and I blame Robert.Woodcock@csiro.au for this.
186         """
187         def __init__(self,domain,D=None,Y=None,q=None,r=None):
188             """
189             initialize the problem
190    
191             @param domain: domain of the PDE.
192             @type domain: L{Domain}
193             @param D: coefficient of the solution.
194             @type D: C{float}, C{int}, L{NumArray}, L{Data}
195             @param Y: right hand side
196             @type Y: C{float}, C{int}, L{NumArray}, L{Data}
197             @param q: location of constraints
198             @type q: C{float}, C{int}, L{NumArray}, L{Data}
199             @param r: value of solution at locations of constraints
200             @type r: C{float}, C{int}, L{NumArray}, L{Data}
201             """
202             self.__domain=domain
203             self.__D=D
204             self.__Y=Y
205             self.__q=q
206             self.__r=r
207             self.__u=None
208             self.__function_space=escript.Solution(self.__domain)
209         def setReducedOn(self):
210             """
211             sets the L{FunctionSpace} of the solution to L{ReducedSolution}
212             """
213             self.__function_space=escript.ReducedSolution(self.__domain)
214             self.__u=None
215    
216         def setReducedOff(self):
217             """
218             sets the L{FunctionSpace} of the solution to L{Solution}
219             """
220             self.__function_space=escript.Solution(self.__domain)
221             self.__u=None
222            
223         def setValue(self,D=None,Y=None,q=None,r=None):
224             """
225             assigns values to the parameters.
226    
227             @param D: coefficient of the solution.
228             @type D: C{float}, C{int}, L{NumArray}, L{Data}
229             @param Y: right hand side
230             @type Y: C{float}, C{int}, L{NumArray}, L{Data}
231             @param q: location of constraints
232             @type q: C{float}, C{int}, L{NumArray}, L{Data}
233             @param r: value of solution at locations of constraints
234             @type r: C{float}, C{int}, L{NumArray}, L{Data}
235             """
236             if not D==None:
237                self.__D=D
238                self.__u=None
239             if not Y==None:
240                self.__Y=Y
241                self.__u=None
242             if not q==None:
243                self.__q=q
244                self.__u=None
245             if not r==None:
246                self.__r=r
247                self.__u=None
248    
249         def getSolution(self):
250             """
251             returns the solution
252            
253             @return: the solution of the problem
254             @rtype: L{Data} object in the L{FunctionSpace} L{Solution} or L{ReducedSolution}.
255             """
256             if self.__u==None:
257                if self.__D==None:
258                   raise ValueError,"coefficient D is undefined"
259                D=escript.Data(self.__D,self.__function_space)
260                if D.getRank()>1:
261                   raise ValueError,"coefficient D must have rank 0 or 1"
262                if self.__Y==None:
263                   self.__u=escript.Data(0.,D.getShape(),self.__function_space)
264                else:
265                   self.__u=util.quotient(self.__Y,D)
266                if not self.__q==None:
267                    q=util.wherePositive(escript.Data(self.__q,self.__function_space))
268                    self.__u*=(1.-q)
269                    if not self.__r==None: self.__u+=q*self.__r
270             return self.__u
271                
272  class Locator:  class Locator:
273       """       """
274       Locator provides access to the values of data objects at a given       Locator provides access to the values of data objects at a given
# Line 91  class Locator: Line 289  class Locator:
289            self.__function_space=where            self.__function_space=where
290         else:         else:
291            self.__function_space=escript.ContinuousFunction(where)            self.__function_space=escript.ContinuousFunction(where)
292         self.__id=util.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()
293    
294       def __str__(self):       def __str__(self):
295         """         """

Legend:
Removed from v.149  
changed lines
  Added in v.609

  ViewVC Help
Powered by ViewVC 1.1.26