escript/py_src/pdetools.py

# $Id$

"""
Provides some tools related to PDEs. 

Currently includes:
    - Projector - to project a discontinuous
"""

import escript
import linearPDEs
import numarray
import util

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):
    """
    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 = linearPDEs.LinearPDE(domain)
    if fast:
      self.__pde.setSolverMethod(linearPDEs.LinearPDE.LUMPING)
    self.__pde.setSymmetryOn()
    self.__pde.setReducedOrderTo(reduce)
    self.__pde.setValue(D = 1.)
    return

  def __del__(self):
    return

  def __call__(self, input_data):
    """
    Projects input_data onto a continuous function

    @param input_data: The input_data to be projected.
    """
    out=escript.Data(0.,input_data.getShape(),what=escript.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 access to 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,escript.FunctionSpace):
          self.__function_space=where
       else:
          self.__function_space=escript.ContinuousFunction(where)
       self.__id=util.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,escript.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

# vim: expandtab shiftwidth=4:
1	jgs	121	# $Id$
2
3			"""
4	jgs	149	Provides some tools related to PDEs.
5	jgs	121
6	jgs	149	Currently includes:
7			- Projector - to project a discontinuous
8	jgs	121	"""
9
10	jgs	149	import escript
11			import linearPDEs
12	jgs	121	import numarray
13	jgs	149	import util
14	jgs	121
15			class Projector:
16	jgs	149	"""
17			The Projector is a factory which projects a discontiuous function onto a
18			continuous function on the a given domain.
19			"""
20	jgs	121	def __init__(self, domain, reduce = True, fast=True):
21			"""
22	jgs	149	Create a continuous function space projector for a domain.
23	jgs	121
24	jgs	149	@param domain: Domain of the projection.
25			@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)
27	jgs	121	"""
28	jgs	149	self.__pde = linearPDEs.LinearPDE(domain)
29	jgs	148	if fast:
30	jgs	149	self.__pde.setSolverMethod(linearPDEs.LinearPDE.LUMPING)
31	jgs	121	self.__pde.setSymmetryOn()
32			self.__pde.setReducedOrderTo(reduce)
33			self.__pde.setValue(D = 1.)
34			return
35
36			def __del__(self):
37			return
38
39			def __call__(self, input_data):
40			"""
41	jgs	149	Projects input_data onto a continuous function
42	jgs	121
43	jgs	149	@param input_data: The input_data to be projected.
44	jgs	121	"""
45	jgs	149	out=escript.Data(0.,input_data.getShape(),what=escript.ContinuousFunction(self.__pde.getDomain()))
46	jgs	121	if input_data.getRank()==0:
47			self.__pde.setValue(Y = input_data)
48			out=self.__pde.getSolution()
49			elif input_data.getRank()==1:
50			for i0 in range(input_data.getShape()[0]):
51			self.__pde.setValue(Y = input_data[i0])
52			out[i0]=self.__pde.getSolution()
53			elif input_data.getRank()==2:
54			for i0 in range(input_data.getShape()[0]):
55			for i1 in range(input_data.getShape()[1]):
56			self.__pde.setValue(Y = input_data[i0,i1])
57			out[i0,i1]=self.__pde.getSolution()
58			elif input_data.getRank()==3:
59			for i0 in range(input_data.getShape()[0]):
60			for i1 in range(input_data.getShape()[1]):
61			for i2 in range(input_data.getShape()[2]):
62			self.__pde.setValue(Y = input_data[i0,i1,i2])
63			out[i0,i1,i2]=self.__pde.getSolution()
64			else:
65			for i0 in range(input_data.getShape()[0]):
66			for i1 in range(input_data.getShape()[1]):
67			for i2 in range(input_data.getShape()[2]):
68			for i3 in range(input_data.getShape()[3]):
69			self.__pde.setValue(Y = input_data[i0,i1,i2,i3])
70			out[i0,i1,i2,i3]=self.__pde.getSolution()
71			return out
72
73
74	jgs	147	class Locator:
75			"""
76	jgs	149	Locator provides access to the values of data objects at a given
77			spatial coordinate x.
78
79			In fact, a Locator object finds the sample in the set of samples of a
80			given function space or domain where which is closest to the given
81			point x.
82	jgs	147	"""
83
84			def __init__(self,where,x=numarray.zeros((3,))):
85	jgs	149	"""
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	jgs	147	self.__function_space=where
92	jgs	121	else:
93	jgs	149	self.__function_space=escript.ContinuousFunction(where)
94			self.__id=util.length(x[:self.__function_space.getDim()]-self.__function_space.getX()).mindp()
95	jgs	121
96	jgs	147	def __str__(self):
97	jgs	149	"""
98			Returns the coordinates of the Locator as a string.
99			"""
100	jgs	147	return "<Locator %s>"%str(self.getX())
101	jgs	121
102	jgs	147	def getFunctionSpace(self):
103	jgs	149	"""
104			Returns the function space of the Locator.
105			"""
106	jgs	147	return self.__function_space
107
108	jgs	121	def getId(self):
109	jgs	149	"""
110			Returns the identifier of the location.
111			"""
112	jgs	121	return self.__id
113
114	jgs	147	def getX(self):
115	jgs	149	"""
116			Returns the exact coordinates of the Locator.
117			"""
118	jgs	147	return self(self.getFunctionSpace().getX())
119	jgs	121
120	jgs	147	def __call__(self,data):
121	jgs	149	"""
122			Returns the value of data at the Locator of a Data object otherwise
123			the object is returned.
124			"""
125	jgs	147	return self.getValue(data)
126	jgs	121
127	jgs	147	def getValue(self,data):
128	jgs	149	"""
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	jgs	147	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	jgs	149
144			# vim: expandtab shiftwidth=4:
Name	Value
svn:eol-style	native
svn:keywords	Author Date Id Revision