/[escript]/trunk/escript/py_src/pdetools.py
ViewVC logotype

Annotation of /trunk/escript/py_src/pdetools.py

Parent Directory Parent Directory | Revision Log Revision Log


Revision 121 - (hide annotations)
Fri May 6 04:26:16 2005 UTC (14 years, 2 months ago) by jgs
Original Path: trunk/esys2/escript/py_src/pdetools.py
File MIME type: text/x-python
File size: 5259 byte(s)
Merge of development branch back to main trunk on 2005-05-06

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     from esys.escript import *
12     from esys.linearPDEs import LinearPDE
13     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     self.__pde.setLumping(fast)
27     self.__pde.setSymmetryOn()
28     self.__pde.setReducedOrderTo(reduce)
29     self.__pde.setValue(D = 1.)
30     return
31    
32     def __del__(self):
33     return
34    
35     def __call__(self, input_data):
36     """
37     @brief projects input_data onto a continuous function
38    
39     @param input_data The input_data to be projected.
40     """
41     out=Data(0.,input_data.getShape(),what=ContinuousFunction(self.__pde.getDomain()))
42     if input_data.getRank()==0:
43     self.__pde.setValue(Y = input_data)
44     out=self.__pde.getSolution()
45     elif input_data.getRank()==1:
46     for i0 in range(input_data.getShape()[0]):
47     self.__pde.setValue(Y = input_data[i0])
48     out[i0]=self.__pde.getSolution()
49     elif input_data.getRank()==2:
50     for i0 in range(input_data.getShape()[0]):
51     for i1 in range(input_data.getShape()[1]):
52     self.__pde.setValue(Y = input_data[i0,i1])
53     out[i0,i1]=self.__pde.getSolution()
54     elif input_data.getRank()==3:
55     for i0 in range(input_data.getShape()[0]):
56     for i1 in range(input_data.getShape()[1]):
57     for i2 in range(input_data.getShape()[2]):
58     self.__pde.setValue(Y = input_data[i0,i1,i2])
59     out[i0,i1,i2]=self.__pde.getSolution()
60     else:
61     for i0 in range(input_data.getShape()[0]):
62     for i1 in range(input_data.getShape()[1]):
63     for i2 in range(input_data.getShape()[2]):
64     for i3 in range(input_data.getShape()[3]):
65     self.__pde.setValue(Y = input_data[i0,i1,i2,i3])
66     out[i0,i1,i2,i3]=self.__pde.getSolution()
67     return out
68    
69    
70     class Location:
71     """Location provides a factory to access the values of data objects at a given spatial coordinate x.
72     In fact, a Location object finds the sample in the set of samples of a given function space or domain which is closest
73     to the given point x. """
74    
75     def __init__(self,x,where):
76     """initializes a Location to access values in Data objects on the Doamin or FunctionSpace where for the sample point which
77     closest to the given point x"""
78     if isinstance(where,FunctionSpace):
79     self.__where=where
80     else:
81     self.__where=ContinuousFunction(where)
82     self.__id=length(x-self.__where.getX()).minarg()
83    
84     def getValue(self,data):
85     """returns the value of data at the Location at a numarray object"""
86     if isinstance(data,Data):
87     return data
88     else:
89     if not data.getFunctionSpace()==self.getFunctionSpace():
90     raise ValueError,"function space of data obejct does not match function space of Location"
91     else:
92     return data.getValue(self.getId())
93     def getX(self):
94     """returns the exact coordinates of the Location"""
95     return self.getValue(self.getFunctionSpace().getX())
96    
97     def getId(self):
98     """returns the identifier of the location"""
99     return self.__id
100    
101     def getFunctionSpace(self):
102     """returns the function space of the Location"""
103     return self.__function_space
104    
105     def __str__(self):
106     """returns the coordinates of the Location as a string"""
107     return str(self.getX())
108    
109     def testProjector(domain):
110     """runs a few test of the Projector factory and returns the largest error plus a description of the test this error occured"""
111     error_max=0.
112     error_text=""
113     x=ContinuousFunction(domain).getX()
114     for f in [True,False]:
115     p=Projector(domain,reduce=False,fast=f)
116     for r in range(5):
117     text="range %s , fast=%s"%(r,f)
118     if r==0:
119     td_ref=x[0]
120     elif r==1:
121     td_ref=x
122     elif r==2:
123     td_ref=[[11.,12.],[21,22.]]*(x[0]+x[1])
124     elif r==3:
125     td_ref=[[[111.,112.],[121,122.]],[[211.,212.],[221,222.]]]*(x[0]+x[1])
126     elif r==3:
127     td_ref=[[[[1111.,1112.],[1121,1122.]],[[1211.,1212.],[1221,1222.]]],[[[2111.,2112.],[2121,2122.]],[[2211.,2212.],[2221,2222.]]]]*(x[0]+x[1])
128     td=p(td_ref.interpolate(Function(domain)))
129     er=Lsup(td-td_ref)/Lsup(td_ref)
130     print text," error = ",er
131     if er>error_max:
132     error_max=er
133     error_text=text
134     return error_max,error_text
135    
136    
137     # this should be removed later
138     if __name__=="__main__":
139     from esys.finley import Rectangle
140     txt=testProjector(Rectangle(56,61))
141     print "test Projector: ",txt
142    
143    
144    

Properties

Name Value
svn:eol-style native
svn:keywords Author Date Id Revision

  ViewVC Help
Powered by ViewVC 1.1.26