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 getValue(self):
      out=self.__v_mem[0]
      if len(out)==1:
          return out[0]
      else:
          return out

  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	def getValue(self):
53	out=self.__v_mem[0]
54	if len(out)==1:
55	return out[0]
56	else:
57	return out
58
59	def checkin(self,dt,*values):
60	"""
61	adds new values to the manager. the p+1 last value get lost
62	"""
63	o=min(self.__order+1,self.__p)
64	self.__order=min(self.__order+1,self.__p)
65	v_mem_new=[values]
66	dt_mem_new=[dt]
67	for i in range(o-1):
68	v_mem_new.append(self.__v_mem[i])
69	dt_mem_new.append(self.__dt_mem[i])
70	v_mem_new.append(self.__v_mem[o-1])
71	self.__order=o
72	self.__v_mem=v_mem_new
73	self.__dt_mem=dt_mem_new
74	self.__t+=dt
75
76	def extrapolate(self,dt):
77	"""
78	extrapolates to dt forward in time.
79	"""
80	if self.__order==0:
81	out=self.__v_mem[0]
82	else:
83	out=[]
84	for i in range(self.__num_val):
85	out.append((1.+dt/self.__dt_mem[0])self.__v_mem[0][i]-dt/self.__dt_mem[0]self.__v_mem[1][i])
86
87	if len(out)==0:
88	return None
89	elif len(out)==1:
90	return out[0]
91	else:
92	return out
93
94
95	class Projector:
96	"""
97	The Projector is a factory which projects a discontiuous function onto a
98	continuous function on the a given domain.
99	"""
100	def __init__(self, domain, reduce = True, fast=True):
101	"""
102	Create a continuous function space projector for a domain.
103
104	@param domain: Domain of the projection.
105	@param reduce: Flag to reduce projection order (default is True)
106	@param fast: Flag to use a fast method based on matrix lumping (default is true)
107	"""
108	self.__pde = linearPDEs.LinearPDE(domain)
109	if fast:
110	self.__pde.setSolverMethod(linearPDEs.LinearPDE.LUMPING)
111	self.__pde.setSymmetryOn()
112	self.__pde.setReducedOrderTo(reduce)
113	self.__pde.setValue(D = 1.)
114	return
115
116	def __del__(self):
117	return
118
119	def __call__(self, input_data):
120	"""
121	Projects input_data onto a continuous function
122
123	@param input_data: The input_data to be projected.
124	"""
125	out=escript.Data(0.,input_data.getShape(),what=escript.ContinuousFunction(self.__pde.getDomain()))
126	if input_data.getRank()==0:
127	self.__pde.setValue(Y = input_data)
128	out=self.__pde.getSolution()
129	elif input_data.getRank()==1:
130	for i0 in range(input_data.getShape()[0]):
131	self.__pde.setValue(Y = input_data[i0])
132	out[i0]=self.__pde.getSolution()
133	elif input_data.getRank()==2:
134	for i0 in range(input_data.getShape()[0]):
135	for i1 in range(input_data.getShape()[1]):
136	self.__pde.setValue(Y = input_data[i0,i1])
137	out[i0,i1]=self.__pde.getSolution()
138	elif input_data.getRank()==3:
139	for i0 in range(input_data.getShape()[0]):
140	for i1 in range(input_data.getShape()[1]):
141	for i2 in range(input_data.getShape()[2]):
142	self.__pde.setValue(Y = input_data[i0,i1,i2])
143	out[i0,i1,i2]=self.__pde.getSolution()
144	else:
145	for i0 in range(input_data.getShape()[0]):
146	for i1 in range(input_data.getShape()[1]):
147	for i2 in range(input_data.getShape()[2]):
148	for i3 in range(input_data.getShape()[3]):
149	self.__pde.setValue(Y = input_data[i0,i1,i2,i3])
150	out[i0,i1,i2,i3]=self.__pde.getSolution()
151	return out
152
153
154	class Locator:
155	"""
156	Locator provides access to the values of data objects at a given
157	spatial coordinate x.
158
159	In fact, a Locator object finds the sample in the set of samples of a
160	given function space or domain where which is closest to the given
161	point x.
162	"""
163
164	def __init__(self,where,x=numarray.zeros((3,))):
165	"""
166	Initializes a Locator to access values in Data objects on the Doamin
167	or FunctionSpace where for the sample point which
168	closest to the given point x.
169	"""
170	if isinstance(where,escript.FunctionSpace):
171	self.__function_space=where
172	else:
173	self.__function_space=escript.ContinuousFunction(where)
174	self.__id=util.length(x[:self.__function_space.getDim()]-self.__function_space.getX()).mindp()
175
176	def __str__(self):
177	"""
178	Returns the coordinates of the Locator as a string.
179	"""
180	return "<Locator %s>"%str(self.getX())
181
182	def getFunctionSpace(self):
183	"""
184	Returns the function space of the Locator.
185	"""
186	return self.__function_space
187
188	def getId(self):
189	"""
190	Returns the identifier of the location.
191	"""
192	return self.__id
193
194	def getX(self):
195	"""
196	Returns the exact coordinates of the Locator.
197	"""
198	return self(self.getFunctionSpace().getX())
199
200	def __call__(self,data):
201	"""
202	Returns the value of data at the Locator of a Data object otherwise
203	the object is returned.
204	"""
205	return self.getValue(data)
206
207	def getValue(self,data):
208	"""
209	Returns the value of data at the Locator if data is a Data object
210	otherwise the object is returned.
211	"""
212	if isinstance(data,escript.Data):
213	if data.getFunctionSpace()==self.getFunctionSpace():
214	out=data.convertToNumArrayFromDPNo(self.getId()[0],self.getId()[1])
215	else:
216	out=data.interpolate(self.getFunctionSpace()).convertToNumArrayFromDPNo(self.getId()[0],self.getId()[1])
217	if data.getRank()==0:
218	return out[0]
219	else:
220	return out
221	else:
222	return data
223
224	# vim: expandtab shiftwidth=4:
Name	Value
svn:eol-style	native
svn:keywords	Author Date Id Revision