/[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 121 by jgs, Fri May 6 04:26:16 2005 UTC trunk/escript/py_src/pdetools.py revision 351 by gross, Tue Dec 13 09:12:15 2005 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.
   
      Projector - to project a discontinuous  
   
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  from esys.escript import *  import escript
13  from esys.linearPDEs import LinearPDE  import linearPDEs
14  import numarray  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    
53      def checkin(self,dt,*values):
54          """
55          adds new values to the manager. the p+1 last value get lost
56          """
57          o=min(self.__order+1,self.__p)
58          self.__order=min(self.__order+1,self.__p)
59          v_mem_new=[values]
60          dt_mem_new=[dt]
61          for i in range(o-1):
62             v_mem_new.append(self.__v_mem[i])
63             dt_mem_new.append(self.__dt_mem[i])
64          v_mem_new.append(self.__v_mem[o-1])
65          self.__order=o
66          self.__v_mem=v_mem_new
67          self.__dt_mem=dt_mem_new
68          self.__t+=dt
69    
70      def extrapolate(self,dt):
71          """
72          extrapolates to dt forward in time.
73          """
74          if self.__order==0:
75             out=self.__v_mem[0]
76          else:
77            out=[]
78            for i in range(self.__num_val):
79               out.append((1.+dt/self.__dt_mem[0])*self.__v_mem[0][i]-dt/self.__dt_mem[0]*self.__v_mem[1][i])
80    
81          if len(out)==0:
82             return None
83          elif len(out)==1:
84             return out[0]
85          else:
86             return out
87    
88  class Projector:  class Projector:
89    """The Projector is a factory which projects a discontiuous function onto a continuous function on the a given domain"""    """
90      The Projector is a factory which projects a discontiuous function onto a
91      continuous function on the a given domain.
92      """
93    def __init__(self, domain, reduce = True, fast=True):    def __init__(self, domain, reduce = True, fast=True):
94      """      """
95      @brief Create a continuous function space projector for a domain.      Create a continuous function space projector for a domain.
96    
97      @param domain Domain of the projection.      @param domain: Domain of the projection.
98      @param reduce Flag to reduce projection order (default is True)      @param reduce: Flag to reduce projection order (default is True)
99      @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)
100      """      """
101      self.__pde = LinearPDE(domain)      self.__pde = linearPDEs.LinearPDE(domain)
102      self.__pde.setLumping(fast)      if fast:
103          self.__pde.setSolverMethod(linearPDEs.LinearPDE.LUMPING)
104      self.__pde.setSymmetryOn()      self.__pde.setSymmetryOn()
105      self.__pde.setReducedOrderTo(reduce)      self.__pde.setReducedOrderTo(reduce)
106      self.__pde.setValue(D = 1.)      self.__pde.setValue(D = 1.)
# Line 34  class Projector: Line 111  class Projector:
111    
112    def __call__(self, input_data):    def __call__(self, input_data):
113      """      """
114      @brief projects input_data onto a continuous function      Projects input_data onto a continuous function
115    
116      @param input_data  The input_data to be projected.      @param input_data: The input_data to be projected.
117      """      """
118      out=Data(0.,input_data.getShape(),what=ContinuousFunction(self.__pde.getDomain()))      out=escript.Data(0.,input_data.getShape(),what=escript.ContinuousFunction(self.__pde.getDomain()))
119      if input_data.getRank()==0:      if input_data.getRank()==0:
120          self.__pde.setValue(Y = input_data)          self.__pde.setValue(Y = input_data)
121          out=self.__pde.getSolution()          out=self.__pde.getSolution()
# Line 67  class Projector: Line 144  class Projector:
144      return out      return out
145    
146    
147  class Location:  class Locator:
148       """Location provides a factory to access the values of data objects at a given spatial coordinate x.       """
149          In fact, a Location object finds the sample in the set of samples of a given function space or domain which is closest       Locator provides access to the values of data objects at a given
150          to the given point x. """       spatial coordinate x.  
151        
152       def __init__(self,x,where):       In fact, a Locator object finds the sample in the set of samples of a
153         """initializes a Location to access values in Data objects on the Doamin or FunctionSpace where for the sample point which       given function space or domain where which is closest to the given
154            closest to the given point x"""       point x.
155         if isinstance(where,FunctionSpace):       """
156            self.__where=where  
157         def __init__(self,where,x=numarray.zeros((3,))):
158           """
159           Initializes a Locator to access values in Data objects on the Doamin
160           or FunctionSpace where for the sample point which
161           closest to the given point x.
162           """
163           if isinstance(where,escript.FunctionSpace):
164              self.__function_space=where
165         else:         else:
166            self.__where=ContinuousFunction(where)            self.__function_space=escript.ContinuousFunction(where)
167         self.__id=length(x-self.__where.getX()).minarg()         self.__id=util.length(x[:self.__function_space.getDim()]-self.__function_space.getX()).mindp()
168    
169       def getValue(self,data):       def __str__(self):
170          """returns the value of data at the Location at a numarray object"""         """
171          if isinstance(data,Data):         Returns the coordinates of the Locator as a string.
172             return data         """
173          else:         return "<Locator %s>"%str(self.getX())
            if not data.getFunctionSpace()==self.getFunctionSpace():  
              raise ValueError,"function space of data obejct does not match function space of Location"  
            else:  
              return data.getValue(self.getId())  
      def getX(self):  
         """returns the exact coordinates of the Location"""  
         return self.getValue(self.getFunctionSpace().getX())  
   
      def getId(self):  
         """returns the identifier of the location"""  
         return self.__id  
174    
175       def getFunctionSpace(self):       def getFunctionSpace(self):
176          """returns the function space of the Location"""          """
177        Returns the function space of the Locator.
178        """
179          return self.__function_space          return self.__function_space
180    
181       def __str__(self):       def getId(self):
182         """returns the coordinates of the Location as a string"""          """
183         return str(self.getX())      Returns the identifier of the location.
184        """
185  def testProjector(domain):          return self.__id
       """runs a few test of the Projector factory and returns the largest error plus a description of the test this error occured"""  
       error_max=0.  
       error_text=""  
       x=ContinuousFunction(domain).getX()  
       for f in [True,False]:  
          p=Projector(domain,reduce=False,fast=f)  
          for r in range(5):  
             text="range %s , fast=%s"%(r,f)  
             if r==0:  
                td_ref=x[0]  
             elif r==1:  
                td_ref=x  
             elif r==2:  
                td_ref=[[11.,12.],[21,22.]]*(x[0]+x[1])  
             elif r==3:  
                td_ref=[[[111.,112.],[121,122.]],[[211.,212.],[221,222.]]]*(x[0]+x[1])  
             elif r==3:  
                td_ref=[[[[1111.,1112.],[1121,1122.]],[[1211.,1212.],[1221,1222.]]],[[[2111.,2112.],[2121,2122.]],[[2211.,2212.],[2221,2222.]]]]*(x[0]+x[1])  
             td=p(td_ref.interpolate(Function(domain)))  
             er=Lsup(td-td_ref)/Lsup(td_ref)  
             print text," error = ",er  
             if er>error_max:  
                  error_max=er  
                  error_text=text  
       return error_max,error_text  
   
   
 # this should be removed later  
 if __name__=="__main__":  
     from esys.finley import Rectangle  
     txt=testProjector(Rectangle(56,61))  
     print "test Projector: ",txt  
186    
187         def getX(self):
188            """
189        Returns the exact coordinates of the Locator.
190        """
191            return self(self.getFunctionSpace().getX())
192    
193         def __call__(self,data):
194            """
195        Returns the value of data at the Locator of a Data object otherwise
196        the object is returned.
197        """
198            return self.getValue(data)
199    
200         def getValue(self,data):
201            """
202        Returns the value of data at the Locator if data is a Data object
203        otherwise the object is returned.
204        """
205            if isinstance(data,escript.Data):
206               if data.getFunctionSpace()==self.getFunctionSpace():
207                 out=data.convertToNumArrayFromDPNo(self.getId()[0],self.getId()[1])
208               else:
209                 out=data.interpolate(self.getFunctionSpace()).convertToNumArrayFromDPNo(self.getId()[0],self.getId()[1])
210               if data.getRank()==0:
211                  return out[0]
212               else:
213                  return out
214            else:
215               return data
216    
217    # vim: expandtab shiftwidth=4:

Legend:
Removed from v.121  
changed lines
  Added in v.351

  ViewVC Help
Powered by ViewVC 1.1.26