escript/py_src/pdetools.py

# $Id$

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

     Projector - to project a discontinuous


"""

from escript import *
from 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 Locator:
     """

        Locator provides a  access the values of data objects at a given spatial coordinate x.
        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. 

     """

     def __init__(self,where,x=numarray.zeros((3,))):
       """initializes a Locator 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.__function_space=where
       else:
          self.__function_space=ContinuousFunction(where)
       self.__id=length(x[:self.__function_space.getDim()]-self.__function_space.getX()).mindp()

     def __str__(self):
       """returns the coordinates of the Locator as a string"""
       return "<Locator %s>"%str(self.getX())

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

     def getId(self):
        """returns the identifier of the location"""
        return self.__id

     def getX(self):
        """returns the exact coordinates of the Locator"""
        return self(self.getFunctionSpace().getX())

     def __call__(self,data):
        """returns the value of data at the Locator of a Data object otherwise the object is returned."""
        return self.getValue(data)

     def getValue(self,data):
        """returns the value of data at the Locator if data is a Data object otherwise the object is returned."""
        if isinstance(data,Data):
           if data.getFunctionSpace()==self.getFunctionSpace():
             out=data.convertToNumArrayFromDPNo(self.getId()[0],self.getId()[1])
           else:
             out=data.interpolate(self.getFunctionSpace()).convertToNumArrayFromDPNo(self.getId()[0],self.getId()[1])
           if data.getRank()==0:
              return out[0]
           else:
              return out
        else:
           return data
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	jgs	146	from escript import *
12			from linearPDEs import LinearPDE
13	jgs	121	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	jgs	147	class Locator:
71			"""
72	jgs	121
73	jgs	147	Locator provides a access the values of data objects at a given spatial coordinate x.
74			In fact, a Locator object finds the sample in the set of samples of a given function space or domain where which is closest
75			to the given point x.
76
77			"""
78
79			def __init__(self,where,x=numarray.zeros((3,))):
80			"""initializes a Locator to access values in Data objects on the Doamin or FunctionSpace where for the sample point which
81	jgs	121	closest to the given point x"""
82			if isinstance(where,FunctionSpace):
83	jgs	147	self.__function_space=where
84	jgs	121	else:
85	jgs	147	self.__function_space=ContinuousFunction(where)
86			self.__id=length(x[:self.__function_space.getDim()]-self.__function_space.getX()).mindp()
87	jgs	121
88	jgs	147	def __str__(self):
89			"""returns the coordinates of the Locator as a string"""
90			return "<Locator %s>"%str(self.getX())
91	jgs	121
92	jgs	147	def getFunctionSpace(self):
93			"""returns the function space of the Locator"""
94			return self.__function_space
95
96	jgs	121	def getId(self):
97			"""returns the identifier of the location"""
98			return self.__id
99
100	jgs	147	def getX(self):
101			"""returns the exact coordinates of the Locator"""
102			return self(self.getFunctionSpace().getX())
103	jgs	121
104	jgs	147	def __call__(self,data):
105			"""returns the value of data at the Locator of a Data object otherwise the object is returned."""
106			return self.getValue(data)
107	jgs	121
108	jgs	147	def getValue(self,data):
109			"""returns the value of data at the Locator if data is a Data object otherwise the object is returned."""
110			if isinstance(data,Data):
111			if data.getFunctionSpace()==self.getFunctionSpace():
112			out=data.convertToNumArrayFromDPNo(self.getId()[0],self.getId()[1])
113			else:
114			out=data.interpolate(self.getFunctionSpace()).convertToNumArrayFromDPNo(self.getId()[0],self.getId()[1])
115			if data.getRank()==0:
116			return out[0]
117			else:
118			return out
119			else:
120			return data
Name	Value
svn:eol-style	native
svn:keywords	Author Date Id Revision