/[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

revision 121 by jgs, Fri May 6 04:26:16 2005 UTC revision 149 by jgs, Thu Sep 1 03:31:39 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  """  """
9    
10  from esys.escript import *  import escript
11  from esys.linearPDEs import LinearPDE  import linearPDEs
12  import numarray  import numarray
13    import util
14    
15  class Projector:  class Projector:
16    """The Projector is a factory which projects a discontiuous function onto a continuous function on the a given domain"""    """
17      The Projector is a factory which projects a discontiuous function onto a
18      continuous function on the a given domain.
19      """
20    def __init__(self, domain, reduce = True, fast=True):    def __init__(self, domain, reduce = True, fast=True):
21      """      """
22      @brief Create a continuous function space projector for a domain.      Create a continuous function space projector for a domain.
23    
24      @param domain Domain of the projection.      @param domain: Domain of the projection.
25      @param reduce Flag to reduce projection order (default is True)      @param reduce: Flag to reduce projection order (default is True)
26      @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)
27      """      """
28      self.__pde = LinearPDE(domain)      self.__pde = linearPDEs.LinearPDE(domain)
29      self.__pde.setLumping(fast)      if fast:
30          self.__pde.setSolverMethod(linearPDEs.LinearPDE.LUMPING)
31      self.__pde.setSymmetryOn()      self.__pde.setSymmetryOn()
32      self.__pde.setReducedOrderTo(reduce)      self.__pde.setReducedOrderTo(reduce)
33      self.__pde.setValue(D = 1.)      self.__pde.setValue(D = 1.)
# Line 34  class Projector: Line 38  class Projector:
38    
39    def __call__(self, input_data):    def __call__(self, input_data):
40      """      """
41      @brief projects input_data onto a continuous function      Projects input_data onto a continuous function
42    
43      @param input_data  The input_data to be projected.      @param input_data: The input_data to be projected.
44      """      """
45      out=Data(0.,input_data.getShape(),what=ContinuousFunction(self.__pde.getDomain()))      out=escript.Data(0.,input_data.getShape(),what=escript.ContinuousFunction(self.__pde.getDomain()))
46      if input_data.getRank()==0:      if input_data.getRank()==0:
47          self.__pde.setValue(Y = input_data)          self.__pde.setValue(Y = input_data)
48          out=self.__pde.getSolution()          out=self.__pde.getSolution()
# Line 67  class Projector: Line 71  class Projector:
71      return out      return out
72    
73    
74  class Location:  class Locator:
75       """Location provides a factory to access the values of data objects at a given spatial coordinate x.       """
76          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
77          to the given point x. """       spatial coordinate x.  
78        
79       def __init__(self,x,where):       In fact, a Locator object finds the sample in the set of samples of a
80         """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
81            closest to the given point x"""       point x.
82         if isinstance(where,FunctionSpace):       """
83            self.__where=where  
84         def __init__(self,where,x=numarray.zeros((3,))):
85           """
86           Initializes a Locator to access values in Data objects on the Doamin
87           or FunctionSpace where for the sample point which
88           closest to the given point x.
89           """
90           if isinstance(where,escript.FunctionSpace):
91              self.__function_space=where
92         else:         else:
93            self.__where=ContinuousFunction(where)            self.__function_space=escript.ContinuousFunction(where)
94         self.__id=length(x-self.__where.getX()).minarg()         self.__id=util.length(x[:self.__function_space.getDim()]-self.__function_space.getX()).mindp()
95    
96       def getValue(self,data):       def __str__(self):
97          """returns the value of data at the Location at a numarray object"""         """
98          if isinstance(data,Data):         Returns the coordinates of the Locator as a string.
99             return data         """
100          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  
101    
102       def getFunctionSpace(self):       def getFunctionSpace(self):
103          """returns the function space of the Location"""          """
104        Returns the function space of the Locator.
105        """
106          return self.__function_space          return self.__function_space
107    
108       def __str__(self):       def getId(self):
109         """returns the coordinates of the Location as a string"""          """
110         return str(self.getX())      Returns the identifier of the location.
111        """
112  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  
113    
114         def getX(self):
115            """
116        Returns the exact coordinates of the Locator.
117        """
118            return self(self.getFunctionSpace().getX())
119    
120         def __call__(self,data):
121            """
122        Returns the value of data at the Locator of a Data object otherwise
123        the object is returned.
124        """
125            return self.getValue(data)
126    
127         def getValue(self,data):
128            """
129        Returns the value of data at the Locator if data is a Data object
130        otherwise the object is returned.
131        """
132            if isinstance(data,escript.Data):
133               if data.getFunctionSpace()==self.getFunctionSpace():
134                 out=data.convertToNumArrayFromDPNo(self.getId()[0],self.getId()[1])
135               else:
136                 out=data.interpolate(self.getFunctionSpace()).convertToNumArrayFromDPNo(self.getId()[0],self.getId()[1])
137               if data.getRank()==0:
138                  return out[0]
139               else:
140                  return out
141            else:
142               return data
143    
144    # vim: expandtab shiftwidth=4:

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

  ViewVC Help
Powered by ViewVC 1.1.26