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	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	jgs	148	if fast:
27			self.__pde.setSolverMethod(LinearPDE.LUMPING)
28	jgs	121	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	jgs	147	class Locator:
72			"""
73	jgs	121
74	jgs	147	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	jgs	121	closest to the given point x"""
83			if isinstance(where,FunctionSpace):
84	jgs	147	self.__function_space=where
85	jgs	121	else:
86	jgs	147	self.__function_space=ContinuousFunction(where)
87			self.__id=length(x[:self.__function_space.getDim()]-self.__function_space.getX()).mindp()
88	jgs	121
89	jgs	147	def __str__(self):
90			"""returns the coordinates of the Locator as a string"""
91			return "<Locator %s>"%str(self.getX())
92	jgs	121
93	jgs	147	def getFunctionSpace(self):
94			"""returns the function space of the Locator"""
95			return self.__function_space
96
97	jgs	121	def getId(self):
98			"""returns the identifier of the location"""
99			return self.__id
100
101	jgs	147	def getX(self):
102			"""returns the exact coordinates of the Locator"""
103			return self(self.getFunctionSpace().getX())
104	jgs	121
105	jgs	147	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	jgs	121
109	jgs	147	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