escript/py_src/pdetools.py

# $Id$

"""
provides a some tools related to PDEs currently includes:

     Projector - to project a discontinuous


"""

from esys.escript import *
from esys.linearPDEs import LinearPDE
import numarray

class Projector:
  """The Projector is a factory which projects a discontiuous function onto a continuous function on the a given domain"""
  def __init__(self, domain, reduce = True, fast=True):
    """
    @brief Create a continuous function space projector for a domain.

    @param domain Domain of the projection.
    @param reduce Flag to reduce projection order (default is True)
    @param fast Flag to use a fast method based on matrix lumping (default is true)
    """
    self.__pde = LinearPDE(domain)
    self.__pde.setLumping(fast)
    self.__pde.setSymmetryOn()
    self.__pde.setReducedOrderTo(reduce)
    self.__pde.setValue(D = 1.)
    return

  def __del__(self):
    return

  def __call__(self, input_data):
    """
    @brief projects input_data onto a continuous function

    @param input_data  The input_data to be projected.
    """
    out=Data(0.,input_data.getShape(),what=ContinuousFunction(self.__pde.getDomain()))
    if input_data.getRank()==0:
        self.__pde.setValue(Y = input_data)
        out=self.__pde.getSolution()
    elif input_data.getRank()==1:
        for i0 in range(input_data.getShape()[0]):
           self.__pde.setValue(Y = input_data[i0])
           out[i0]=self.__pde.getSolution()
    elif input_data.getRank()==2:
        for i0 in range(input_data.getShape()[0]):
           for i1 in range(input_data.getShape()[1]):
              self.__pde.setValue(Y = input_data[i0,i1])
              out[i0,i1]=self.__pde.getSolution()
    elif input_data.getRank()==3:
        for i0 in range(input_data.getShape()[0]):
           for i1 in range(input_data.getShape()[1]):
              for i2 in range(input_data.getShape()[2]):
                 self.__pde.setValue(Y = input_data[i0,i1,i2])
                 out[i0,i1,i2]=self.__pde.getSolution()
    else:
        for i0 in range(input_data.getShape()[0]):
           for i1 in range(input_data.getShape()[1]):
              for i2 in range(input_data.getShape()[2]):
                 for i3 in range(input_data.getShape()[3]):
                    self.__pde.setValue(Y = input_data[i0,i1,i2,i3])
                    out[i0,i1,i2,i3]=self.__pde.getSolution()
    return out


class Location:
     """Location provides a factory to access the values of data objects at a given spatial coordinate x.
        In fact, a Location object finds the sample in the set of samples of a given function space or domain which is closest
        to the given point x. """

     def __init__(self,x,where):
       """initializes a Location to access values in Data objects on the Doamin or FunctionSpace where for the sample point which
          closest to the given point x"""
       if isinstance(where,FunctionSpace):
          self.__where=where
       else:
          self.__where=ContinuousFunction(where)
       self.__id=length(x-self.__where.getX()).minarg()

     def getValue(self,data):
        """returns the value of data at the Location at a numarray object"""
        if isinstance(data,Data):
           return data
        else:
           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

     def getFunctionSpace(self):
        """returns the function space of the Location"""
        return self.__function_space

     def __str__(self):
       """returns the coordinates of the Location as a string"""
       return str(self.getX())

def testProjector(domain):
      """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


1	jgs	121	# $Id$
2
3			"""
4			provides a some tools related to PDEs currently includes:
5
6			Projector - to project a discontinuous
7
8
9			"""
10
11			from esys.escript import *
12			from esys.linearPDEs import LinearPDE
13			import numarray
14
15			class Projector:
16			"""The Projector is a factory which projects a discontiuous function onto a continuous function on the a given domain"""
17			def __init__(self, domain, reduce = True, fast=True):
18			"""
19			@brief Create a continuous function space projector for a domain.
20
21			@param domain Domain of the projection.
22			@param reduce Flag to reduce projection order (default is True)
23			@param fast Flag to use a fast method based on matrix lumping (default is true)
24			"""
25			self.__pde = LinearPDE(domain)
26			self.__pde.setLumping(fast)
27			self.__pde.setSymmetryOn()
28			self.__pde.setReducedOrderTo(reduce)
29			self.__pde.setValue(D = 1.)
30			return
31
32			def __del__(self):
33			return
34
35			def __call__(self, input_data):
36			"""
37			@brief projects input_data onto a continuous function
38
39			@param input_data The input_data to be projected.
40			"""
41			out=Data(0.,input_data.getShape(),what=ContinuousFunction(self.__pde.getDomain()))
42			if input_data.getRank()==0:
43			self.__pde.setValue(Y = input_data)
44			out=self.__pde.getSolution()
45			elif input_data.getRank()==1:
46			for i0 in range(input_data.getShape()[0]):
47			self.__pde.setValue(Y = input_data[i0])
48			out[i0]=self.__pde.getSolution()
49			elif input_data.getRank()==2:
50			for i0 in range(input_data.getShape()[0]):
51			for i1 in range(input_data.getShape()[1]):
52			self.__pde.setValue(Y = input_data[i0,i1])
53			out[i0,i1]=self.__pde.getSolution()
54			elif input_data.getRank()==3:
55			for i0 in range(input_data.getShape()[0]):
56			for i1 in range(input_data.getShape()[1]):
57			for i2 in range(input_data.getShape()[2]):
58			self.__pde.setValue(Y = input_data[i0,i1,i2])
59			out[i0,i1,i2]=self.__pde.getSolution()
60			else:
61			for i0 in range(input_data.getShape()[0]):
62			for i1 in range(input_data.getShape()[1]):
63			for i2 in range(input_data.getShape()[2]):
64			for i3 in range(input_data.getShape()[3]):
65			self.__pde.setValue(Y = input_data[i0,i1,i2,i3])
66			out[i0,i1,i2,i3]=self.__pde.getSolution()
67			return out
68
69
70			class Location:
71			"""Location provides a factory to access the values of data objects at a given spatial coordinate x.
72			In fact, a Location object finds the sample in the set of samples of a given function space or domain which is closest
73			to the given point x. """
74
75			def __init__(self,x,where):
76			"""initializes a Location to access values in Data objects on the Doamin or FunctionSpace where for the sample point which
77			closest to the given point x"""
78			if isinstance(where,FunctionSpace):
79			self.__where=where
80			else:
81			self.__where=ContinuousFunction(where)
82			self.__id=length(x-self.__where.getX()).minarg()
83
84			def getValue(self,data):
85			"""returns the value of data at the Location at a numarray object"""
86			if isinstance(data,Data):
87			return data
88			else:
89			if not data.getFunctionSpace()==self.getFunctionSpace():
90			raise ValueError,"function space of data obejct does not match function space of Location"
91			else:
92			return data.getValue(self.getId())
93			def getX(self):
94			"""returns the exact coordinates of the Location"""
95			return self.getValue(self.getFunctionSpace().getX())
96
97			def getId(self):
98			"""returns the identifier of the location"""
99			return self.__id
100
101			def getFunctionSpace(self):
102			"""returns the function space of the Location"""
103			return self.__function_space
104
105			def __str__(self):
106			"""returns the coordinates of the Location as a string"""
107			return str(self.getX())
108
109			def testProjector(domain):
110			"""runs a few test of the Projector factory and returns the largest error plus a description of the test this error occured"""
111			error_max=0.
112			error_text=""
113			x=ContinuousFunction(domain).getX()
114			for f in [True,False]:
115			p=Projector(domain,reduce=False,fast=f)
116			for r in range(5):
117			text="range %s , fast=%s"%(r,f)
118			if r==0:
119			td_ref=x[0]
120			elif r==1:
121			td_ref=x
122			elif r==2:
123			td_ref=[[11.,12.],[21,22.]]*(x[0]+x[1])
124			elif r==3:
125			td_ref=[[[111.,112.],[121,122.]],[[211.,212.],[221,222.]]]*(x[0]+x[1])
126			elif r==3:
127			td_ref=[[[[1111.,1112.],[1121,1122.]],[[1211.,1212.],[1221,1222.]]],[[[2111.,2112.],[2121,2122.]],[[2211.,2212.],[2221,2222.]]]]*(x[0]+x[1])
128			td=p(td_ref.interpolate(Function(domain)))
129			er=Lsup(td-td_ref)/Lsup(td_ref)
130			print text," error = ",er
131			if er>error_max:
132			error_max=er
133			error_text=text
134			return error_max,error_text
135
136
137			# this should be removed later
138			if __name__=="__main__":
139			from esys.finley import Rectangle
140			txt=testProjector(Rectangle(56,61))
141			print "test Projector: ",txt
142
143
144
Name	Value
svn:eol-style	native
svn:keywords	Author Date Id Revision