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)
    if fast:
      self.__pde.setSolverMethod(LinearPDE.LUMPING)
    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	# $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 escript import *
12	from 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	if fast:
27	self.__pde.setSolverMethod(LinearPDE.LUMPING)
28	self.__pde.setSymmetryOn()
29	self.__pde.setReducedOrderTo(reduce)
30	self.__pde.setValue(D = 1.)
31	return
32
33	def __del__(self):
34	return
35
36	def __call__(self, input_data):
37	"""
38	@brief projects input_data onto a continuous function
39
40	@param input_data The input_data to be projected.
41	"""
42	out=Data(0.,input_data.getShape(),what=ContinuousFunction(self.__pde.getDomain()))
43	if input_data.getRank()==0:
44	self.__pde.setValue(Y = input_data)
45	out=self.__pde.getSolution()
46	elif input_data.getRank()==1:
47	for i0 in range(input_data.getShape()[0]):
48	self.__pde.setValue(Y = input_data[i0])
49	out[i0]=self.__pde.getSolution()
50	elif input_data.getRank()==2:
51	for i0 in range(input_data.getShape()[0]):
52	for i1 in range(input_data.getShape()[1]):
53	self.__pde.setValue(Y = input_data[i0,i1])
54	out[i0,i1]=self.__pde.getSolution()
55	elif input_data.getRank()==3:
56	for i0 in range(input_data.getShape()[0]):
57	for i1 in range(input_data.getShape()[1]):
58	for i2 in range(input_data.getShape()[2]):
59	self.__pde.setValue(Y = input_data[i0,i1,i2])
60	out[i0,i1,i2]=self.__pde.getSolution()
61	else:
62	for i0 in range(input_data.getShape()[0]):
63	for i1 in range(input_data.getShape()[1]):
64	for i2 in range(input_data.getShape()[2]):
65	for i3 in range(input_data.getShape()[3]):
66	self.__pde.setValue(Y = input_data[i0,i1,i2,i3])
67	out[i0,i1,i2,i3]=self.__pde.getSolution()
68	return out
69
70
71	class Locator:
72	"""
73
74	Locator provides a access the values of data objects at a given spatial coordinate x.
75	In fact, a Locator object finds the sample in the set of samples of a given function space or domain where which is closest
76	to the given point x.
77
78	"""
79
80	def __init__(self,where,x=numarray.zeros((3,))):
81	"""initializes a Locator to access values in Data objects on the Doamin or FunctionSpace where for the sample point which
82	closest to the given point x"""
83	if isinstance(where,FunctionSpace):
84	self.__function_space=where
85	else:
86	self.__function_space=ContinuousFunction(where)
87	self.__id=length(x[:self.__function_space.getDim()]-self.__function_space.getX()).mindp()
88
89	def __str__(self):
90	"""returns the coordinates of the Locator as a string"""
91	return "<Locator %s>"%str(self.getX())
92
93	def getFunctionSpace(self):
94	"""returns the function space of the Locator"""
95	return self.__function_space
96
97	def getId(self):
98	"""returns the identifier of the location"""
99	return self.__id
100
101	def getX(self):
102	"""returns the exact coordinates of the Locator"""
103	return self(self.getFunctionSpace().getX())
104
105	def __call__(self,data):
106	"""returns the value of data at the Locator of a Data object otherwise the object is returned."""
107	return self.getValue(data)
108
109	def getValue(self,data):
110	"""returns the value of data at the Locator if data is a Data object otherwise the object is returned."""
111	if isinstance(data,Data):
112	if data.getFunctionSpace()==self.getFunctionSpace():
113	out=data.convertToNumArrayFromDPNo(self.getId()[0],self.getId()[1])
114	else:
115	out=data.interpolate(self.getFunctionSpace()).convertToNumArrayFromDPNo(self.getId()[0],self.getId()[1])
116	if data.getRank()==0:
117	return out[0]
118	else:
119	return out
120	else:
121	return data
Name	Value
svn:eol-style	native
svn:keywords	Author Date Id Revision