/[escript]/trunk/escript/test/python/test_linearPDEs.py
ViewVC logotype

Diff of /trunk/escript/test/python/test_linearPDEs.py

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

revision 2323 by gross, Thu Mar 19 04:23:32 2009 UTC revision 2325 by gross, Thu Mar 19 05:32:13 2009 UTC
# Line 44  The tests must be linked with a Domain c Line 44  The tests must be linked with a Domain c
44    
45  __author__="Lutz Gross, l.gross@uq.edu.au"  __author__="Lutz Gross, l.gross@uq.edu.au"
46    
47  from esys.escript.util import Lsup,kronecker,interpolate,whereZero  from esys.escript.util import Lsup,kronecker,interpolate,whereZero, outer, swap_axes
48  from esys.escript import Function,FunctionOnBoundary,FunctionOnContactZero,Solution,ReducedSolution,Vector,ContinuousFunction,Scalar, ReducedFunction,ReducedFunctionOnBoundary,ReducedFunctionOnContactZero,Data, Tensor4, Tensor  from esys.escript import Function,FunctionOnBoundary,FunctionOnContactZero,Solution,ReducedSolution,Vector,ContinuousFunction,Scalar, ReducedFunction,ReducedFunctionOnBoundary,ReducedFunctionOnContactZero,Data, Tensor4, Tensor
49  from esys.escript.linearPDEs import LinearPDE,IllegalCoefficientValue,Poisson, IllegalCoefficientFunctionSpace, TransportPDE, IllegalCoefficient, Helmholtz  from esys.escript.linearPDEs import LinearPDE,IllegalCoefficientValue,Poisson, IllegalCoefficientFunctionSpace, TransportPDE, IllegalCoefficient, Helmholtz, LameEquation
50  import numarray  import numarray
51  import unittest  import unittest
52    
# Line 62  class Test_linearPDEs(unittest.TestCase) Line 62  class Test_linearPDEs(unittest.TestCase)
62          if tol==None: tol=self.TOL          if tol==None: tol=self.TOL
63          return Lsup(arg-ref_arg)<=tol*Lsup(ref_arg)          return Lsup(arg-ref_arg)<=tol*Lsup(ref_arg)
64            
65    class Test_LameEquation(Test_linearPDEs):
66    
67        def test_config(self):
68            mypde=LameEquation(self.domain,debug=self.DEBUG)
69            d=self.domain.getDim()
70            self.failUnlessEqual((mypde.getNumEquations(),mypde.getNumSolutions(),mypde.isSymmetric()),(d,d,True),"set up incorrect")
71    
72        def test_setCoefficient_q(self):
73            mypde=LameEquation(self.domain,debug=self.DEBUG)
74            x=self.domain.getX()
75            mypde.setValue(q=x)
76    
77            q_ref=interpolate(x,Solution(self.domain))
78            self.failUnless(self.check(mypde.getCoefficient("A"),0),"A is not 0")
79            self.failUnless(mypde.getCoefficient("B").isEmpty(),"B is not empty")
80            self.failUnless(mypde.getCoefficient("C").isEmpty(),"C is not empty")
81            self.failUnless(mypde.getCoefficient("D").isEmpty(),"D is not empty")
82            self.failUnless(mypde.getCoefficient("X").isEmpty(),"X is not empty")
83            self.failUnless(mypde.getCoefficient("Y").isEmpty(),"Y is not empty")
84            self.failUnless(mypde.getCoefficient("y").isEmpty(),"y is not empty")
85            self.failUnless(mypde.getCoefficient("d").isEmpty(),"d is not empty")
86            self.failUnless(mypde.getCoefficient("d_contact").isEmpty(),"d_contact is not empty")
87            self.failUnless(mypde.getCoefficient("y_contact").isEmpty(),"y_contact is not empty")
88            self.failUnless(mypde.getCoefficient("A_reduced").isEmpty(),"A_reduced is not empty")
89            self.failUnless(mypde.getCoefficient("B_reduced").isEmpty(),"B_reduced is not empty")
90            self.failUnless(mypde.getCoefficient("C_reduced").isEmpty(),"C_reduced is not empty")
91            self.failUnless(mypde.getCoefficient("D_reduced").isEmpty(),"D_reduced is not empty")
92            self.failUnless(mypde.getCoefficient("X_reduced").isEmpty(),"X_reduced is not empty")
93            self.failUnless(mypde.getCoefficient("Y_reduced").isEmpty(),"Y_reduced is not empty")
94            self.failUnless(mypde.getCoefficient("y_reduced").isEmpty(),"y_reduced is not empty")
95            self.failUnless(mypde.getCoefficient("d_reduced").isEmpty(),"d_reduced is not empty")
96            self.failUnless(mypde.getCoefficient("d_contact_reduced").isEmpty(),"d_contact_reduced is not empty")
97            self.failUnless(mypde.getCoefficient("y_contact_reduced").isEmpty(),"y_contact_reduced is not empty")
98            self.failUnless(self.check(mypde.getCoefficient("q"),q_ref),"q is not empty")
99            self.failUnless(mypde.getCoefficient("r").isEmpty(),"r is not empty")
100    
101        def test_setCoefficient_r(self):
102            mypde=LameEquation(self.domain,debug=self.DEBUG)
103            x=self.domain.getX()
104            mypde.setValue(r=x)
105    
106            r_ref=interpolate(x,Solution(self.domain))
107            self.failUnless(self.check(mypde.getCoefficient("A"),0),"A is not 0")
108            self.failUnless(mypde.getCoefficient("B").isEmpty(),"B is not empty")
109            self.failUnless(mypde.getCoefficient("C").isEmpty(),"C is not empty")
110            self.failUnless(mypde.getCoefficient("D").isEmpty(),"D is not empty")
111            self.failUnless(mypde.getCoefficient("X").isEmpty(),"X is not empty")
112            self.failUnless(mypde.getCoefficient("Y").isEmpty(),"Y is not empty")
113            self.failUnless(mypde.getCoefficient("y").isEmpty(),"y is not empty")
114            self.failUnless(mypde.getCoefficient("d").isEmpty(),"d is not empty")
115            self.failUnless(mypde.getCoefficient("d_contact").isEmpty(),"d_contact is not empty")
116            self.failUnless(mypde.getCoefficient("y_contact").isEmpty(),"y_contact is not empty")
117            self.failUnless(mypde.getCoefficient("A_reduced").isEmpty(),"A_reduced is not empty")
118            self.failUnless(mypde.getCoefficient("B_reduced").isEmpty(),"B_reduced is not empty")
119            self.failUnless(mypde.getCoefficient("C_reduced").isEmpty(),"C_reduced is not empty")
120            self.failUnless(mypde.getCoefficient("D_reduced").isEmpty(),"D_reduced is not empty")
121            self.failUnless(mypde.getCoefficient("X_reduced").isEmpty(),"X_reduced is not empty")
122            self.failUnless(mypde.getCoefficient("Y_reduced").isEmpty(),"Y_reduced is not empty")
123            self.failUnless(mypde.getCoefficient("y_reduced").isEmpty(),"y_reduced is not empty")
124            self.failUnless(mypde.getCoefficient("d_reduced").isEmpty(),"d_reduced is not empty")
125            self.failUnless(mypde.getCoefficient("d_contact_reduced").isEmpty(),"d_contact_reduced is not empty")
126            self.failUnless(mypde.getCoefficient("y_contact_reduced").isEmpty(),"y_contact_reduced is not empty")
127            self.failUnless(self.check(mypde.getCoefficient("r"),r_ref),"r is nor x")
128            self.failUnless(mypde.getCoefficient("q").isEmpty(),"q is not empty")
129    
130    
131        def test_setCoefficient_F(self):
132            mypde=LameEquation(self.domain,debug=self.DEBUG)
133            x=self.domain.getX()
134            mypde.setValue(F=x)
135    
136            Y_ref=interpolate(x,Function(self.domain))
137            self.failUnless(self.check(mypde.getCoefficient("A"),0),"A is not 0")
138            self.failUnless(mypde.getCoefficient("B").isEmpty(),"B is not empty")
139            self.failUnless(mypde.getCoefficient("C").isEmpty(),"C is not empty")
140            self.failUnless(mypde.getCoefficient("D").isEmpty(),"D is not empty")
141            self.failUnless(mypde.getCoefficient("X").isEmpty(),"X is not empty")
142            self.failUnless(self.check(mypde.getCoefficient("Y"),Y_ref),"Y is not x")
143            self.failUnless(mypde.getCoefficient("y").isEmpty(),"y is not empty")
144            self.failUnless(mypde.getCoefficient("d").isEmpty(),"d is not empty")
145            self.failUnless(mypde.getCoefficient("d_contact").isEmpty(),"d_contact is not empty")
146            self.failUnless(mypde.getCoefficient("y_contact").isEmpty(),"y_contact is not empty")
147            self.failUnless(mypde.getCoefficient("A_reduced").isEmpty(),"A_reduced is not empty")
148            self.failUnless(mypde.getCoefficient("B_reduced").isEmpty(),"B_reduced is not empty")
149            self.failUnless(mypde.getCoefficient("C_reduced").isEmpty(),"C_reduced is not empty")
150            self.failUnless(mypde.getCoefficient("D_reduced").isEmpty(),"D_reduced is not empty")
151            self.failUnless(mypde.getCoefficient("X_reduced").isEmpty(),"X_reduced is not empty")
152            self.failUnless(mypde.getCoefficient("Y_reduced").isEmpty(),"Y_reduced is not empty")
153            self.failUnless(mypde.getCoefficient("y_reduced").isEmpty(),"y_reduced is not empty")
154            self.failUnless(mypde.getCoefficient("d_reduced").isEmpty(),"d_reduced is not empty")
155            self.failUnless(mypde.getCoefficient("d_contact_reduced").isEmpty(),"d_contact_reduced is not empty")
156            self.failUnless(mypde.getCoefficient("y_contact_reduced").isEmpty(),"y_contact_reduced is not empty")
157            self.failUnless(mypde.getCoefficient("q").isEmpty(),"q is not empty")
158            self.failUnless(mypde.getCoefficient("r").isEmpty(),"r is not empty")
159    
160        def test_setCoefficient_f(self):
161            mypde=LameEquation(self.domain,debug=self.DEBUG)
162            x=self.domain.getX()
163            mypde.setValue(f=x)
164    
165            y_ref=interpolate(x,FunctionOnBoundary(self.domain))
166            self.failUnless(self.check(mypde.getCoefficient("A"),0),"A is not 0")
167            self.failUnless(mypde.getCoefficient("B").isEmpty(),"B is not empty")
168            self.failUnless(mypde.getCoefficient("C").isEmpty(),"C is not empty")
169            self.failUnless(mypde.getCoefficient("D").isEmpty(),"D is not empty")
170            self.failUnless(mypde.getCoefficient("X").isEmpty(),"X is not empty")
171            self.failUnless(mypde.getCoefficient("Y").isEmpty(),"Y is not empty")
172            self.failUnless(mypde.getCoefficient("d").isEmpty(),"d is not empty")
173            self.failUnless(self.check(mypde.getCoefficient("y"),y_ref),"d is not x[0]")
174            self.failUnless(mypde.getCoefficient("d_contact").isEmpty(),"d_contact is not empty")
175            self.failUnless(mypde.getCoefficient("y_contact").isEmpty(),"y_contact is not empty")
176            self.failUnless(mypde.getCoefficient("A_reduced").isEmpty(),"A_reduced is not empty")
177            self.failUnless(mypde.getCoefficient("B_reduced").isEmpty(),"B_reduced is not empty")
178            self.failUnless(mypde.getCoefficient("C_reduced").isEmpty(),"C_reduced is not empty")
179            self.failUnless(mypde.getCoefficient("D_reduced").isEmpty(),"D_reduced is not empty")
180            self.failUnless(mypde.getCoefficient("X_reduced").isEmpty(),"X_reduced is not empty")
181            self.failUnless(mypde.getCoefficient("Y_reduced").isEmpty(),"X_reduced is not empty")
182            self.failUnless(mypde.getCoefficient("y_reduced").isEmpty(),"y_reduced is not empty")
183            self.failUnless(mypde.getCoefficient("d_reduced").isEmpty(),"d_reduced is not empty")
184            self.failUnless(mypde.getCoefficient("d_contact_reduced").isEmpty(),"d_contact_reduced is not empty")
185            self.failUnless(mypde.getCoefficient("y_contact_reduced").isEmpty(),"y_contact_reduced is not empty")
186            self.failUnless(mypde.getCoefficient("q").isEmpty(),"q is not empty")
187            self.failUnless(mypde.getCoefficient("r").isEmpty(),"r is not empty")
188    
189        def test_setCoefficient_sigma(self):
190            mypde=LameEquation(self.domain,debug=self.DEBUG)
191            x=self.domain.getX()
192            mypde.setValue(sigma=outer(x,x))
193    
194            X_ref=interpolate(outer(x,x),Function(self.domain))
195            self.failUnless(self.check(mypde.getCoefficient("A"),0),"A is not 0")
196            self.failUnless(mypde.getCoefficient("B").isEmpty(),"B is not empty")
197            self.failUnless(mypde.getCoefficient("C").isEmpty(),"C is not empty")
198            self.failUnless(mypde.getCoefficient("D").isEmpty(),"D is not empty")
199            self.failUnless(self.check(mypde.getCoefficient("X"),X_ref),"X is not x X x")
200            self.failUnless(mypde.getCoefficient("Y").isEmpty(),"Y is not empty")
201            self.failUnless(mypde.getCoefficient("d").isEmpty(),"d is not empty")
202            self.failUnless(mypde.getCoefficient("y").isEmpty(),"y is not empty")
203            self.failUnless(mypde.getCoefficient("d_contact").isEmpty(),"d_contact is not empty")
204            self.failUnless(mypde.getCoefficient("y_contact").isEmpty(),"y_contact is not empty")
205            self.failUnless(mypde.getCoefficient("A_reduced").isEmpty(),"A_reduced is not empty")
206            self.failUnless(mypde.getCoefficient("B_reduced").isEmpty(),"B_reduced is not empty")
207            self.failUnless(mypde.getCoefficient("C_reduced").isEmpty(),"C_reduced is not empty")
208            self.failUnless(mypde.getCoefficient("D_reduced").isEmpty(),"D_reduced is not empty")
209            self.failUnless(mypde.getCoefficient("X_reduced").isEmpty(),"X_reduced is not empty")
210            self.failUnless(mypde.getCoefficient("Y_reduced").isEmpty(),"X_reduced is not empty")
211            self.failUnless(mypde.getCoefficient("y_reduced").isEmpty(),"y_reduced is not empty")
212            self.failUnless(mypde.getCoefficient("d_reduced").isEmpty(),"d_reduced is not empty")
213            self.failUnless(mypde.getCoefficient("d_contact_reduced").isEmpty(),"d_contact_reduced is not empty")
214            self.failUnless(mypde.getCoefficient("y_contact_reduced").isEmpty(),"y_contact_reduced is not empty")
215            self.failUnless(mypde.getCoefficient("q").isEmpty(),"q is not empty")
216            self.failUnless(mypde.getCoefficient("r").isEmpty(),"r is not empty")
217    
218        def test_setCoefficient_lambda(self):
219            mypde=LameEquation(self.domain,debug=self.DEBUG)
220            x=self.domain.getX()
221            mypde.setValue(lame_lambda=x[0])
222    
223    
224            k3=kronecker(Function(self.domain))
225            k3Xk3=outer(k3,k3)
226            A_ref=x[0]*k3Xk3
227    
228            self.failUnless(self.check(mypde.getCoefficient("A"),A_ref),"A is not kronecker")
229            self.failUnless(mypde.getCoefficient("B").isEmpty(),"B is not empty")
230            self.failUnless(mypde.getCoefficient("C").isEmpty(),"C is not empty")
231            self.failUnless(mypde.getCoefficient("D").isEmpty(),"D is not empty")
232            self.failUnless(mypde.getCoefficient("X").isEmpty(),"X is not empty")
233            self.failUnless(mypde.getCoefficient("Y").isEmpty(),"Y is not empty")
234            self.failUnless(mypde.getCoefficient("d").isEmpty(),"d is not empty")
235            self.failUnless(mypde.getCoefficient("y").isEmpty(),"y is not empty")
236            self.failUnless(mypde.getCoefficient("d_contact").isEmpty(),"d_contact is not empty")
237            self.failUnless(mypde.getCoefficient("y_contact").isEmpty(),"y_contact is not empty")
238            self.failUnless(mypde.getCoefficient("A_reduced").isEmpty(),"A_reduced is not empty")
239            self.failUnless(mypde.getCoefficient("B_reduced").isEmpty(),"B_reduced is not empty")
240            self.failUnless(mypde.getCoefficient("C_reduced").isEmpty(),"C_reduced is not empty")
241            self.failUnless(mypde.getCoefficient("D_reduced").isEmpty(),"D_reduced is not empty")
242            self.failUnless(mypde.getCoefficient("X_reduced").isEmpty(),"X_reduced is not empty")
243            self.failUnless(mypde.getCoefficient("Y_reduced").isEmpty(),"Y_reduced is not empty")
244            self.failUnless(mypde.getCoefficient("y_reduced").isEmpty(),"y_reduced is not empty")
245            self.failUnless(mypde.getCoefficient("d_reduced").isEmpty(),"d_reduced is not empty")
246            self.failUnless(mypde.getCoefficient("d_contact_reduced").isEmpty(),"d_contact_reduced is not empty")
247            self.failUnless(mypde.getCoefficient("y_contact_reduced").isEmpty(),"y_contact_reduced is not empty")
248            self.failUnless(mypde.getCoefficient("q").isEmpty(),"q is not empty")
249            self.failUnless(mypde.getCoefficient("r").isEmpty(),"r is not empty")
250    
251        def test_setCoefficient_mu(self):
252            mypde=LameEquation(self.domain,debug=self.DEBUG)
253            x=self.domain.getX()
254            mypde.setValue(lame_mu=x[0])
255    
256    
257            k3=kronecker(Function(self.domain))
258            k3Xk3=outer(k3,k3)
259            A_ref=x[0]*(swap_axes(k3Xk3,0,3)+swap_axes(k3Xk3,1,3))
260    
261            self.failUnless(self.check(mypde.getCoefficient("A"),A_ref),"A is not kronecker")
262            self.failUnless(mypde.getCoefficient("B").isEmpty(),"B is not empty")
263            self.failUnless(mypde.getCoefficient("C").isEmpty(),"C is not empty")
264            self.failUnless(mypde.getCoefficient("D").isEmpty(),"D is not empty")
265            self.failUnless(mypde.getCoefficient("X").isEmpty(),"X is not empty")
266            self.failUnless(mypde.getCoefficient("Y").isEmpty(),"Y is not empty")
267            self.failUnless(mypde.getCoefficient("d").isEmpty(),"d is not empty")
268            self.failUnless(mypde.getCoefficient("y").isEmpty(),"y is not empty")
269            self.failUnless(mypde.getCoefficient("d_contact").isEmpty(),"d_contact is not empty")
270            self.failUnless(mypde.getCoefficient("y_contact").isEmpty(),"y_contact is not empty")
271            self.failUnless(mypde.getCoefficient("A_reduced").isEmpty(),"A_reduced is not empty")
272            self.failUnless(mypde.getCoefficient("B_reduced").isEmpty(),"B_reduced is not empty")
273            self.failUnless(mypde.getCoefficient("C_reduced").isEmpty(),"C_reduced is not empty")
274            self.failUnless(mypde.getCoefficient("D_reduced").isEmpty(),"D_reduced is not empty")
275            self.failUnless(mypde.getCoefficient("X_reduced").isEmpty(),"X_reduced is not empty")
276            self.failUnless(mypde.getCoefficient("Y_reduced").isEmpty(),"Y_reduced is not empty")
277            self.failUnless(mypde.getCoefficient("y_reduced").isEmpty(),"y_reduced is not empty")
278            self.failUnless(mypde.getCoefficient("d_reduced").isEmpty(),"d_reduced is not empty")
279            self.failUnless(mypde.getCoefficient("d_contact_reduced").isEmpty(),"d_contact_reduced is not empty")
280            self.failUnless(mypde.getCoefficient("y_contact_reduced").isEmpty(),"y_contact_reduced is not empty")
281            self.failUnless(mypde.getCoefficient("q").isEmpty(),"q is not empty")
282            self.failUnless(mypde.getCoefficient("r").isEmpty(),"r is not empty")
283    
284        def test_setCoefficient_lambdamu(self):
285            mypde=LameEquation(self.domain,debug=self.DEBUG)
286            x=self.domain.getX()
287            mypde.setValue(lame_lambda=x[0], lame_mu=x[1])
288    
289            k3=kronecker(Function(self.domain))
290            k3Xk3=outer(k3,k3)
291            A_ref=x[0]*k3Xk3+x[1]*(swap_axes(k3Xk3,0,3)+swap_axes(k3Xk3,1,3))
292    
293            self.failUnless(self.check(mypde.getCoefficient("A"),A_ref),"A is not kronecker")
294            self.failUnless(mypde.getCoefficient("B").isEmpty(),"B is not empty")
295            self.failUnless(mypde.getCoefficient("C").isEmpty(),"C is not empty")
296            self.failUnless(mypde.getCoefficient("D").isEmpty(),"D is not empty")
297            self.failUnless(mypde.getCoefficient("X").isEmpty(),"X is not empty")
298            self.failUnless(mypde.getCoefficient("Y").isEmpty(),"Y is not empty")
299            self.failUnless(mypde.getCoefficient("d").isEmpty(),"d is not empty")
300            self.failUnless(mypde.getCoefficient("y").isEmpty(),"y is not empty")
301            self.failUnless(mypde.getCoefficient("d_contact").isEmpty(),"d_contact is not empty")
302            self.failUnless(mypde.getCoefficient("y_contact").isEmpty(),"y_contact is not empty")
303            self.failUnless(mypde.getCoefficient("A_reduced").isEmpty(),"A_reduced is not empty")
304            self.failUnless(mypde.getCoefficient("B_reduced").isEmpty(),"B_reduced is not empty")
305            self.failUnless(mypde.getCoefficient("C_reduced").isEmpty(),"C_reduced is not empty")
306            self.failUnless(mypde.getCoefficient("D_reduced").isEmpty(),"D_reduced is not empty")
307            self.failUnless(mypde.getCoefficient("X_reduced").isEmpty(),"X_reduced is not empty")
308            self.failUnless(mypde.getCoefficient("Y_reduced").isEmpty(),"Y_reduced is not empty")
309            self.failUnless(mypde.getCoefficient("y_reduced").isEmpty(),"y_reduced is not empty")
310            self.failUnless(mypde.getCoefficient("d_reduced").isEmpty(),"d_reduced is not empty")
311            self.failUnless(mypde.getCoefficient("d_contact_reduced").isEmpty(),"d_contact_reduced is not empty")
312            self.failUnless(mypde.getCoefficient("y_contact_reduced").isEmpty(),"y_contact_reduced is not empty")
313            self.failUnless(mypde.getCoefficient("q").isEmpty(),"q is not empty")
314            self.failUnless(mypde.getCoefficient("r").isEmpty(),"r is not empty")
315    
316        def test_solve(self):
317           d=self.domain.getDim()
318           mypde=LameEquation(self.domain)
319           cf=ContinuousFunction(self.domain)
320           x=cf.getX()
321           u_ex=x
322           msk=Vector(0.,cf)
323           for i in range(d): msk[i]=whereZero(x[i])
324           mypde.setValue(q=msk,r=u_ex,lame_mu=3,lame_lambda=50,f=(2*3+50*d)*FunctionOnBoundary(self.domain).getNormal())
325    
326           u=mypde.getSolution()
327           self.failUnless(self.check(u,u_ex,10*self.TOL),"incorrect solution")
328    
329  class Test_Helmholtz(Test_linearPDEs):  class Test_Helmholtz(Test_linearPDEs):
330    
331      def test_config(self):      def test_config(self):

Legend:
Removed from v.2323  
changed lines
  Added in v.2325

  ViewVC Help
Powered by ViewVC 1.1.26