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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 146 - (show annotations)
Fri Jul 29 01:44:29 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: 5249 byte(s)
First merge of development branch dev-02 back to main trunk on 2005-07-29

1 # $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 escript import *
12 from 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