escript/py_src/pdetools.py

# $Id$

"""
Provides some tools related to PDEs. 

Currently includes:
    - Projector - to project a discontinuous
    - Locator - to trace values in data objects at a certain location
    - TimeIntegrationManager - to handel extraplotion in time
"""

import escript
import linearPDEs
import numarray
import util

class TimeIntegrationManager:
  """
  a simple mechanism to manage time dependend values. 

  typical usage is:

  dt=0.1 # time increment
  tm=TimeIntegrationManager(inital_value,p=1)
  while t<1.
      v_guess=tm.extrapolate(dt) # extrapolate to t+dt
      v=...
      tm.checkin(dt,v)
      t+=dt

  @remark: currently only p=1 is supported.
  """
  def __init__(self,*inital_values,**kwargs):
     """
     sets up the value manager where inital_value is the initial value and p is order used for extrapolation
     """
     if kwargs.has_key("p"):
            self.__p=kwargs["p"]
     else:
            self.__p=1
     if kwargs.has_key("time"):
            self.__t=kwargs["time"]
     else:
            self.__t=0.
     self.__v_mem=[inital_values]
     self.__order=0
     self.__dt_mem=[]
     self.__num_val=len(inital_values)

  def getTime(self):
      return self.__t

  def checkin(self,dt,*values):
      """
      adds new values to the manager. the p+1 last value get lost
      """
      o=min(self.__order+1,self.__p)
      self.__order=min(self.__order+1,self.__p)
      v_mem_new=[values]
      dt_mem_new=[dt]
      for i in range(o-1):
         v_mem_new.append(self.__v_mem[i])
         dt_mem_new.append(self.__dt_mem[i])
      v_mem_new.append(self.__v_mem[o-1])
      self.__order=o
      self.__v_mem=v_mem_new
      self.__dt_mem=dt_mem_new
      self.__t+=dt

  def extrapolate(self,dt):
      """
      extrapolates to dt forward in time.
      """
      if self.__order==0:
         out=self.__v_mem[0]
      else:
        out=[]
        for i in range(self.__num_val):
           out.append((1.+dt/self.__dt_mem[0])*self.__v_mem[0][i]-dt/self.__dt_mem[0]*self.__v_mem[1][i])

      if len(out)==0:
         return None
      elif len(out)==1:
         return out[0]
      else:
         return out

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	# $Id$
2
3	"""
4	Provides some tools related to PDEs.
5
6	Currently includes:
7	- Projector - to project a discontinuous
8	- Locator - to trace values in data objects at a certain location
9	- TimeIntegrationManager - to handel extraplotion in time
10	"""
11
12	import escript
13	import linearPDEs
14	import numarray
15	import util
16
17	class TimeIntegrationManager:
18	"""
19	a simple mechanism to manage time dependend values.
20
21	typical usage is:
22
23	dt=0.1 # time increment
24	tm=TimeIntegrationManager(inital_value,p=1)
25	while t<1.
26	v_guess=tm.extrapolate(dt) # extrapolate to t+dt
27	v=...
28	tm.checkin(dt,v)
29	t+=dt
30
31	@remark: currently only p=1 is supported.
32	"""
33	def __init__(self,inital_values,*kwargs):
34	"""
35	sets up the value manager where inital_value is the initial value and p is order used for extrapolation
36	"""
37	if kwargs.has_key("p"):
38	self.__p=kwargs["p"]
39	else:
40	self.__p=1
41	if kwargs.has_key("time"):
42	self.__t=kwargs["time"]
43	else:
44	self.__t=0.
45	self.__v_mem=[inital_values]
46	self.__order=0
47	self.__dt_mem=[]
48	self.__num_val=len(inital_values)
49
50	def getTime(self):
51	return self.__t
52
53	def checkin(self,dt,*values):
54	"""
55	adds new values to the manager. the p+1 last value get lost
56	"""
57	o=min(self.__order+1,self.__p)
58	self.__order=min(self.__order+1,self.__p)
59	v_mem_new=[values]
60	dt_mem_new=[dt]
61	for i in range(o-1):
62	v_mem_new.append(self.__v_mem[i])
63	dt_mem_new.append(self.__dt_mem[i])
64	v_mem_new.append(self.__v_mem[o-1])
65	self.__order=o
66	self.__v_mem=v_mem_new
67	self.__dt_mem=dt_mem_new
68	self.__t+=dt
69
70	def extrapolate(self,dt):
71	"""
72	extrapolates to dt forward in time.
73	"""
74	if self.__order==0:
75	out=self.__v_mem[0]
76	else:
77	out=[]
78	for i in range(self.__num_val):
79	out.append((1.+dt/self.__dt_mem[0])self.__v_mem[0][i]-dt/self.__dt_mem[0]self.__v_mem[1][i])
80
81	if len(out)==0:
82	return None
83	elif len(out)==1:
84	return out[0]
85	else:
86	return out
87
88	class Projector:
89	"""
90	The Projector is a factory which projects a discontiuous function onto a
91	continuous function on the a given domain.
92	"""
93	def __init__(self, domain, reduce = True, fast=True):
94	"""
95	Create a continuous function space projector for a domain.
96
97	@param domain: Domain of the projection.
98	@param reduce: Flag to reduce projection order (default is True)
99	@param fast: Flag to use a fast method based on matrix lumping (default is true)
100	"""
101	self.__pde = linearPDEs.LinearPDE(domain)
102	if fast:
103	self.__pde.setSolverMethod(linearPDEs.LinearPDE.LUMPING)
104	self.__pde.setSymmetryOn()
105	self.__pde.setReducedOrderTo(reduce)
106	self.__pde.setValue(D = 1.)
107	return
108
109	def __del__(self):
110	return
111
112	def __call__(self, input_data):
113	"""
114	Projects input_data onto a continuous function
115
116	@param input_data: The input_data to be projected.
117	"""
118	out=escript.Data(0.,input_data.getShape(),what=escript.ContinuousFunction(self.__pde.getDomain()))
119	if input_data.getRank()==0:
120	self.__pde.setValue(Y = input_data)
121	out=self.__pde.getSolution()
122	elif input_data.getRank()==1:
123	for i0 in range(input_data.getShape()[0]):
124	self.__pde.setValue(Y = input_data[i0])
125	out[i0]=self.__pde.getSolution()
126	elif input_data.getRank()==2:
127	for i0 in range(input_data.getShape()[0]):
128	for i1 in range(input_data.getShape()[1]):
129	self.__pde.setValue(Y = input_data[i0,i1])
130	out[i0,i1]=self.__pde.getSolution()
131	elif input_data.getRank()==3:
132	for i0 in range(input_data.getShape()[0]):
133	for i1 in range(input_data.getShape()[1]):
134	for i2 in range(input_data.getShape()[2]):
135	self.__pde.setValue(Y = input_data[i0,i1,i2])
136	out[i0,i1,i2]=self.__pde.getSolution()
137	else:
138	for i0 in range(input_data.getShape()[0]):
139	for i1 in range(input_data.getShape()[1]):
140	for i2 in range(input_data.getShape()[2]):
141	for i3 in range(input_data.getShape()[3]):
142	self.__pde.setValue(Y = input_data[i0,i1,i2,i3])
143	out[i0,i1,i2,i3]=self.__pde.getSolution()
144	return out
145
146
147	class Locator:
148	"""
149	Locator provides access to the values of data objects at a given
150	spatial coordinate x.
151
152	In fact, a Locator object finds the sample in the set of samples of a
153	given function space or domain where which is closest to the given
154	point x.
155	"""
156
157	def __init__(self,where,x=numarray.zeros((3,))):
158	"""
159	Initializes a Locator to access values in Data objects on the Doamin
160	or FunctionSpace where for the sample point which
161	closest to the given point x.
162	"""
163	if isinstance(where,escript.FunctionSpace):
164	self.__function_space=where
165	else:
166	self.__function_space=escript.ContinuousFunction(where)
167	self.__id=util.length(x[:self.__function_space.getDim()]-self.__function_space.getX()).mindp()
168
169	def __str__(self):
170	"""
171	Returns the coordinates of the Locator as a string.
172	"""
173	return "<Locator %s>"%str(self.getX())
174
175	def getFunctionSpace(self):
176	"""
177	Returns the function space of the Locator.
178	"""
179	return self.__function_space
180
181	def getId(self):
182	"""
183	Returns the identifier of the location.
184	"""
185	return self.__id
186
187	def getX(self):
188	"""
189	Returns the exact coordinates of the Locator.
190	"""
191	return self(self.getFunctionSpace().getX())
192
193	def __call__(self,data):
194	"""
195	Returns the value of data at the Locator of a Data object otherwise
196	the object is returned.
197	"""
198	return self.getValue(data)
199
200	def getValue(self,data):
201	"""
202	Returns the value of data at the Locator if data is a Data object
203	otherwise the object is returned.
204	"""
205	if isinstance(data,escript.Data):
206	if data.getFunctionSpace()==self.getFunctionSpace():
207	out=data.convertToNumArrayFromDPNo(self.getId()[0],self.getId()[1])
208	else:
209	out=data.interpolate(self.getFunctionSpace()).convertToNumArrayFromDPNo(self.getId()[0],self.getId()[1])
210	if data.getRank()==0:
211	return out[0]
212	else:
213	return out
214	else:
215	return data
216
217	# vim: expandtab shiftwidth=4:
Name	Value
svn:eol-style	native
svn:keywords	Author Date Id Revision