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

@var __author__: name of author
@var __copyright__: copyrights
@var __license__: licence agreement
@var __url__: url entry point on documentation
@var __version__: version
@var __date__: date of the version
"""

__author__="Lutz Gross, l.gross@uq.edu.au"
__copyright__="""  Copyright (c) 2006 by ACcESS MNRF
                    http://www.access.edu.au
                Primary Business: Queensland, Australia"""
__license__="""Licensed under the Open Software License version 3.0
             http://www.opensource.org/licenses/osl-3.0.php"""
__url__="http://www.iservo.edu.au/esys"
__version__="$Revision$"
__date__="$Date$"


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

  @note: 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(),self.__pde.getFunctionSpaceForSolution())
    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 NoPDE:
     """
     solves the following problem for u:

     M{kronecker[i,j]*D[j]*u[j]=Y[i]} 

     with constraint

     M{u[j]=r[j]}  where M{q[j]>0}

     where D, Y, r and q are given functions of rank 1.

     In the case of scalars this takes the form

     M{D*u=Y} 

     with constraint

     M{u=r}  where M{q>0}

     where D, Y, r and q are given scalar functions.

     The constraint is overwriting any other condition.

     @note: This class is similar to the L{linearPDEs.LinearPDE} class with A=B=C=X=0 but has the intention
            that all input parameter are given in L{Solution} or L{ReducedSolution}. The whole
            thing is a bit strange and I blame Robert.Woodcock@csiro.au for this.
     """
     def __init__(self,domain,D=None,Y=None,q=None,r=None):
         """
         initialize the problem

         @param domain: domain of the PDE.
         @type domain: L{Domain}
         @param D: coefficient of the solution. 
         @type D: C{float}, C{int}, L{numarray.NumArray}, L{Data}
         @param Y: right hand side
         @type Y: C{float}, C{int}, L{numarray.NumArray}, L{Data}
         @param q: location of constraints
         @type q: C{float}, C{int}, L{numarray.NumArray}, L{Data}
         @param r: value of solution at locations of constraints
         @type r: C{float}, C{int}, L{numarray.NumArray}, L{Data}
         """
         self.__domain=domain
         self.__D=D
         self.__Y=Y
         self.__q=q
         self.__r=r
         self.__u=None
         self.__function_space=escript.Solution(self.__domain)
     def setReducedOn(self):
         """
         sets the L{FunctionSpace} of the solution to L{ReducedSolution}
         """
         self.__function_space=escript.ReducedSolution(self.__domain)
         self.__u=None

     def setReducedOff(self):
         """
         sets the L{FunctionSpace} of the solution to L{Solution}
         """
         self.__function_space=escript.Solution(self.__domain)
         self.__u=None
         
     def setValue(self,D=None,Y=None,q=None,r=None):
         """
         assigns values to the parameters.

         @param D: coefficient of the solution. 
         @type D: C{float}, C{int}, L{numarray.NumArray}, L{Data}
         @param Y: right hand side
         @type Y: C{float}, C{int}, L{numarray.NumArray}, L{Data}
         @param q: location of constraints
         @type q: C{float}, C{int}, L{numarray.NumArray}, L{Data}
         @param r: value of solution at locations of constraints
         @type r: C{float}, C{int}, L{numarray.NumArray}, L{Data}
         """
         if not D==None:
            self.__D=D
            self.__u=None
         if not Y==None:
            self.__Y=Y
            self.__u=None
         if not q==None:
            self.__q=q
            self.__u=None
         if not r==None:
            self.__r=r
            self.__u=None

     def getSolution(self):
         """
         returns the solution
        
         @return: the solution of the problem
         @rtype: L{Data} object in the L{FunctionSpace} L{Solution} or L{ReducedSolution}.
         """
         if self.__u==None:
            if self.__D==None:
               raise ValueError,"coefficient D is undefined"
            D=escript.Data(self.__D,self.__function_space)
            if D.getRank()>1:
               raise ValueError,"coefficient D must have rank 0 or 1"
            if self.__Y==None:
               self.__u=escript.Data(0.,D.getShape(),self.__function_space)
            else:
               self.__u=util.quotient(self.__Y,D)
            if not self.__q==None:
                q=util.wherePositive(escript.Data(self.__q,self.__function_space))
                self.__u*=(1.-q)
                if not self.__r==None: self.__u+=q*self.__r
         return self.__u
             
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(self.__function_space.getX()-x[:self.__function_space.getDim()]).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])
             out=data.convertToNumArrayFromDPNo(self.getId()[0],self.getId()[1],self.getId()[2])
           else:
             #out=data.interpolate(self.getFunctionSpace()).convertToNumArrayFromDPNo(self.getId()[0],self.getId()[1])
             out=data.interpolate(self.getFunctionSpace()).convertToNumArrayFromDPNo(self.getId()[0],self.getId()[1],self.getId()[2])
           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	@var __author__: name of author
12	@var __copyright__: copyrights
13	@var __license__: licence agreement
14	@var __url__: url entry point on documentation
15	@var __version__: version
16	@var __date__: date of the version
17	"""
18
19	__author__="Lutz Gross, l.gross@uq.edu.au"
20	__copyright__=""" Copyright (c) 2006 by ACcESS MNRF
21	http://www.access.edu.au
22	Primary Business: Queensland, Australia"""
23	__license__="""Licensed under the Open Software License version 3.0
24	http://www.opensource.org/licenses/osl-3.0.php"""
25	__url__="http://www.iservo.edu.au/esys"
26	__version__="$Revision$"
27	__date__="$Date$"
28
29
30	import escript
31	import linearPDEs
32	import numarray
33	import util
34
35	class TimeIntegrationManager:
36	"""
37	a simple mechanism to manage time dependend values.
38
39	typical usage is::
40
41	dt=0.1 # time increment
42	tm=TimeIntegrationManager(inital_value,p=1)
43	while t<1.
44	v_guess=tm.extrapolate(dt) # extrapolate to t+dt
45	v=...
46	tm.checkin(dt,v)
47	t+=dt
48
49	@note: currently only p=1 is supported.
50	"""
51	def __init__(self,inital_values,*kwargs):
52	"""
53	sets up the value manager where inital_value is the initial value and p is order used for extrapolation
54	"""
55	if kwargs.has_key("p"):
56	self.__p=kwargs["p"]
57	else:
58	self.__p=1
59	if kwargs.has_key("time"):
60	self.__t=kwargs["time"]
61	else:
62	self.__t=0.
63	self.__v_mem=[inital_values]
64	self.__order=0
65	self.__dt_mem=[]
66	self.__num_val=len(inital_values)
67
68	def getTime(self):
69	return self.__t
70	def getValue(self):
71	out=self.__v_mem[0]
72	if len(out)==1:
73	return out[0]
74	else:
75	return out
76
77	def checkin(self,dt,*values):
78	"""
79	adds new values to the manager. the p+1 last value get lost
80	"""
81	o=min(self.__order+1,self.__p)
82	self.__order=min(self.__order+1,self.__p)
83	v_mem_new=[values]
84	dt_mem_new=[dt]
85	for i in range(o-1):
86	v_mem_new.append(self.__v_mem[i])
87	dt_mem_new.append(self.__dt_mem[i])
88	v_mem_new.append(self.__v_mem[o-1])
89	self.__order=o
90	self.__v_mem=v_mem_new
91	self.__dt_mem=dt_mem_new
92	self.__t+=dt
93
94	def extrapolate(self,dt):
95	"""
96	extrapolates to dt forward in time.
97	"""
98	if self.__order==0:
99	out=self.__v_mem[0]
100	else:
101	out=[]
102	for i in range(self.__num_val):
103	out.append((1.+dt/self.__dt_mem[0])self.__v_mem[0][i]-dt/self.__dt_mem[0]self.__v_mem[1][i])
104
105	if len(out)==0:
106	return None
107	elif len(out)==1:
108	return out[0]
109	else:
110	return out
111
112
113	class Projector:
114	"""
115	The Projector is a factory which projects a discontiuous function onto a
116	continuous function on the a given domain.
117	"""
118	def __init__(self, domain, reduce = True, fast=True):
119	"""
120	Create a continuous function space projector for a domain.
121
122	@param domain: Domain of the projection.
123	@param reduce: Flag to reduce projection order (default is True)
124	@param fast: Flag to use a fast method based on matrix lumping (default is true)
125	"""
126	self.__pde = linearPDEs.LinearPDE(domain)
127	if fast:
128	self.__pde.setSolverMethod(linearPDEs.LinearPDE.LUMPING)
129	self.__pde.setSymmetryOn()
130	self.__pde.setReducedOrderTo(reduce)
131	self.__pde.setValue(D = 1.)
132	return
133
134	def __del__(self):
135	return
136
137	def __call__(self, input_data):
138	"""
139	Projects input_data onto a continuous function
140
141	@param input_data: The input_data to be projected.
142	"""
143	out=escript.Data(0.,input_data.getShape(),self.__pde.getFunctionSpaceForSolution())
144	if input_data.getRank()==0:
145	self.__pde.setValue(Y = input_data)
146	out=self.__pde.getSolution()
147	elif input_data.getRank()==1:
148	for i0 in range(input_data.getShape()[0]):
149	self.__pde.setValue(Y = input_data[i0])
150	out[i0]=self.__pde.getSolution()
151	elif input_data.getRank()==2:
152	for i0 in range(input_data.getShape()[0]):
153	for i1 in range(input_data.getShape()[1]):
154	self.__pde.setValue(Y = input_data[i0,i1])
155	out[i0,i1]=self.__pde.getSolution()
156	elif input_data.getRank()==3:
157	for i0 in range(input_data.getShape()[0]):
158	for i1 in range(input_data.getShape()[1]):
159	for i2 in range(input_data.getShape()[2]):
160	self.__pde.setValue(Y = input_data[i0,i1,i2])
161	out[i0,i1,i2]=self.__pde.getSolution()
162	else:
163	for i0 in range(input_data.getShape()[0]):
164	for i1 in range(input_data.getShape()[1]):
165	for i2 in range(input_data.getShape()[2]):
166	for i3 in range(input_data.getShape()[3]):
167	self.__pde.setValue(Y = input_data[i0,i1,i2,i3])
168	out[i0,i1,i2,i3]=self.__pde.getSolution()
169	return out
170
171	class NoPDE:
172	"""
173	solves the following problem for u:
174
175	M{kronecker[i,j]D[j]u[j]=Y[i]}
176
177	with constraint
178
179	M{u[j]=r[j]} where M{q[j]>0}
180
181	where D, Y, r and q are given functions of rank 1.
182
183	In the case of scalars this takes the form
184
185	M{D*u=Y}
186
187	with constraint
188
189	M{u=r} where M{q>0}
190
191	where D, Y, r and q are given scalar functions.
192
193	The constraint is overwriting any other condition.
194
195	@note: This class is similar to the L{linearPDEs.LinearPDE} class with A=B=C=X=0 but has the intention
196	that all input parameter are given in L{Solution} or L{ReducedSolution}. The whole
197	thing is a bit strange and I blame Robert.Woodcock@csiro.au for this.
198	"""
199	def __init__(self,domain,D=None,Y=None,q=None,r=None):
200	"""
201	initialize the problem
202
203	@param domain: domain of the PDE.
204	@type domain: L{Domain}
205	@param D: coefficient of the solution.
206	@type D: C{float}, C{int}, L{numarray.NumArray}, L{Data}
207	@param Y: right hand side
208	@type Y: C{float}, C{int}, L{numarray.NumArray}, L{Data}
209	@param q: location of constraints
210	@type q: C{float}, C{int}, L{numarray.NumArray}, L{Data}
211	@param r: value of solution at locations of constraints
212	@type r: C{float}, C{int}, L{numarray.NumArray}, L{Data}
213	"""
214	self.__domain=domain
215	self.__D=D
216	self.__Y=Y
217	self.__q=q
218	self.__r=r
219	self.__u=None
220	self.__function_space=escript.Solution(self.__domain)
221	def setReducedOn(self):
222	"""
223	sets the L{FunctionSpace} of the solution to L{ReducedSolution}
224	"""
225	self.__function_space=escript.ReducedSolution(self.__domain)
226	self.__u=None
227
228	def setReducedOff(self):
229	"""
230	sets the L{FunctionSpace} of the solution to L{Solution}
231	"""
232	self.__function_space=escript.Solution(self.__domain)
233	self.__u=None
234
235	def setValue(self,D=None,Y=None,q=None,r=None):
236	"""
237	assigns values to the parameters.
238
239	@param D: coefficient of the solution.
240	@type D: C{float}, C{int}, L{numarray.NumArray}, L{Data}
241	@param Y: right hand side
242	@type Y: C{float}, C{int}, L{numarray.NumArray}, L{Data}
243	@param q: location of constraints
244	@type q: C{float}, C{int}, L{numarray.NumArray}, L{Data}
245	@param r: value of solution at locations of constraints
246	@type r: C{float}, C{int}, L{numarray.NumArray}, L{Data}
247	"""
248	if not D==None:
249	self.__D=D
250	self.__u=None
251	if not Y==None:
252	self.__Y=Y
253	self.__u=None
254	if not q==None:
255	self.__q=q
256	self.__u=None
257	if not r==None:
258	self.__r=r
259	self.__u=None
260
261	def getSolution(self):
262	"""
263	returns the solution
264
265	@return: the solution of the problem
266	@rtype: L{Data} object in the L{FunctionSpace} L{Solution} or L{ReducedSolution}.
267	"""
268	if self.__u==None:
269	if self.__D==None:
270	raise ValueError,"coefficient D is undefined"
271	D=escript.Data(self.__D,self.__function_space)
272	if D.getRank()>1:
273	raise ValueError,"coefficient D must have rank 0 or 1"
274	if self.__Y==None:
275	self.__u=escript.Data(0.,D.getShape(),self.__function_space)
276	else:
277	self.__u=util.quotient(self.__Y,D)
278	if not self.__q==None:
279	q=util.wherePositive(escript.Data(self.__q,self.__function_space))
280	self.__u*=(1.-q)
281	if not self.__r==None: self.__u+=q*self.__r
282	return self.__u
283
284	class Locator:
285	"""
286	Locator provides access to the values of data objects at a given
287	spatial coordinate x.
288
289	In fact, a Locator object finds the sample in the set of samples of a
290	given function space or domain where which is closest to the given
291	point x.
292	"""
293
294	def __init__(self,where,x=numarray.zeros((3,))):
295	"""
296	Initializes a Locator to access values in Data objects on the Doamin
297	or FunctionSpace where for the sample point which
298	closest to the given point x.
299	"""
300	if isinstance(where,escript.FunctionSpace):
301	self.__function_space=where
302	else:
303	self.__function_space=escript.ContinuousFunction(where)
304	self.__id=util.length(self.__function_space.getX()-x[:self.__function_space.getDim()]).mindp()
305
306	def __str__(self):
307	"""
308	Returns the coordinates of the Locator as a string.
309	"""
310	return "<Locator %s>"%str(self.getX())
311
312	def getFunctionSpace(self):
313	"""
314	Returns the function space of the Locator.
315	"""
316	return self.__function_space
317
318	def getId(self):
319	"""
320	Returns the identifier of the location.
321	"""
322	return self.__id
323
324	def getX(self):
325	"""
326	Returns the exact coordinates of the Locator.
327	"""
328	return self(self.getFunctionSpace().getX())
329
330	def __call__(self,data):
331	"""
332	Returns the value of data at the Locator of a Data object otherwise
333	the object is returned.
334	"""
335	return self.getValue(data)
336
337	def getValue(self,data):
338	"""
339	Returns the value of data at the Locator if data is a Data object
340	otherwise the object is returned.
341	"""
342	if isinstance(data,escript.Data):
343	if data.getFunctionSpace()==self.getFunctionSpace():
344	#out=data.convertToNumArrayFromDPNo(self.getId()[0],self.getId()[1])
345	out=data.convertToNumArrayFromDPNo(self.getId()[0],self.getId()[1],self.getId()[2])
346	else:
347	#out=data.interpolate(self.getFunctionSpace()).convertToNumArrayFromDPNo(self.getId()[0],self.getId()[1])
348	out=data.interpolate(self.getFunctionSpace()).convertToNumArrayFromDPNo(self.getId()[0],self.getId()[1],self.getId()[2])
349	if data.getRank()==0:
350	return out[0]
351	else:
352	return out
353	else:
354	return data
355
356	# vim: expandtab shiftwidth=4:
Name	Value
svn:eol-style	native
svn:keywords	Author Date Id Revision