/[escript]/trunk/finley/test/python/gentest
ViewVC logotype

Contents of /trunk/finley/test/python/gentest

Parent Directory Parent Directory | Revision Log Revision Log


Revision 766 - (show annotations)
Fri Jun 30 06:44:14 2006 UTC (13 years, 7 months ago) by gross
File size: 31719 byte(s)
modified version of a test generator (private use)
1 #!/usr/bin/python
2 import numarray
3 import random
4 t_prog={}
5 t_prog["2Do1_strong"]=""
6 t_prog["2Do2_strong"]=""
7 t_prog["3Do1_strong"]=""
8 t_prog["3Do2_strong"]=""
9 t_prog["2Do1_weak"]=""
10 t_prog["2Do2_weak"]=""
11 t_prog["3Do1_weak"]=""
12 t_prog["3Do2_weak"]=""
13 t_prog["2Do1_contact"]=""
14 t_prog["2Do2_contact"]=""
15 t_prog["3Do1_contact"]=""
16 t_prog["3Do2_contact"]=""
17 intend=" "
18 max_order=3
19
20 def makeTitle(d,coeffo,o,s,coef,case,typ,body="",mark="",pdeargs="",add_jump=False):
21 if mark=="":
22 name="test_assemblage_%sD_solO%s_coeffO%s_neq%s_%s_%s_type%s"%(d,o,coeffo,s,coef,case,typ)
23 else:
24 name="test_assemblage_%sD_solO%s_coeffO%s_neq%s_%s_%s_type%s_comp%s"%(d,o,coeffo,s,coef,case,typ,mark)
25 key="%sDo%s_%s"%(d,o,typ.lower())
26 t_prog[key]+=intend+"#"+50*"="+"\n"+intend+"def %s(self):\n"%name
27 t_prog[key]+=2*intend+"x=self.domain.getX()\n"
28 if add_jump:
29 t_prog[key]+=2*intend+"jump=Data(0.,(),ContinuousFunction(self.domain))\n"
30 t_prog[key]+=2*intend+"jump.setTaggedValue(2,1.)\n"
31 t_prog[key]+=body
32 t_prog[key]+=2*intend+"pde=LinearPDE(self.domain)\n"
33 t_prog[key]+=2*intend+"pde.setValue(%s)\n"%pdeargs
34 t_prog[key]+=2*intend+"r=pde.getResidual(u)\n"
35 t_prog[key]+=2*intend+"rhs=pde.getRightHandSide()\n"
36 t_prog[key]+=2*intend+"self.failUnless(Lsup(rhs)>0,\"right hand side is zero\")\n"
37 t_prog[key]+=2*intend+"self.failUnless(Lsup(r)<=self.RES_TOL*Lsup(rhs),\"residual is too big\")\n"
38
39 def getRandom():
40 out=0
41 while out==0:
42 out = int(18*random.random())-9
43 return out
44
45 def makeTestSolution(order,ncomps,dim):
46 if dim == 2:
47 out=numarray.zeros((ncomps,max_order+1,max_order+1))
48 for d1 in range(order+1):
49 for d2 in range(order+1):
50 for k in range(ncomps):
51 out[k,d1,d2]=getRandom()
52 if dim == 3:
53 out=numarray.zeros((ncomps,max_order+1,max_order+1,max_order+1))
54 for d1 in range(order+1):
55 for d2 in range(order+1):
56 for d3 in range(order+1):
57 for k in range(ncomps):
58 out[k,d1,d2,d3]=getRandom()
59 return out
60
61 def makeGradient(func,dim):
62 sh=func.shape
63 if len(sh)==dim+1:
64 sh_out=sh[:1]+(dim,)+sh[1:]
65 out=numarray.zeros(sh_out)
66 if dim==2:
67 for d1 in range(sh[1]-1):
68 for d2 in range(sh[2]-1):
69 for k in range(sh[0]):
70 out[k,0,d1,d2]=(d1+1)*func[k,d1+1,d2]
71 out[k,1,d1,d2]=(d2+1)*func[k,d1,d2+1]
72 else:
73 for d1 in range(sh[1]-1):
74 for d2 in range(sh[2]-1):
75 for d3 in range(sh[3]-1):
76 for k in range(sh[0]):
77 out[k,0,d1,d2,d3]=(d1+1)*func[k,d1+1,d2,d3]
78 out[k,1,d1,d2,d3]=(d2+1)*func[k,d1,d2+1,d3]
79 out[k,2,d1,d2,d3]=(d3+1)*func[k,d1,d2,d3+1]
80 elif len(sh)==dim+2:
81 raise RunTimeError,"sdfadsa"
82 sh_out=sh[:2]+(dim,)+sh[2:]
83 out=numarray.zeros(sh_out)
84 else:
85 raise RunTimeError,"sdfadsa"
86 return out
87
88 def makeDiv(func,dim):
89 sh=func.shape
90 if len(sh)==dim+1:
91 sh_out=(1,)+sh[1:]
92 out=numarray.zeros(sh_out)
93 if dim==2:
94 for d1 in range(sh[1]-1):
95 for d2 in range(sh[2]-1):
96 for k in range(sh[0]):
97 out[0,d1,d2]=(d1+1)*func[0,d1+1,d2]+(d2+1)*func[1,d1,d2+1]
98 else:
99 for d1 in range(sh[1]-1):
100 for d2 in range(sh[2]-1):
101 for d3 in range(sh[3]-1):
102 for k in range(sh[0]):
103 out[0,d1,d2,d3]=(d1+1)*func[0,d1+1,d2,d3]+(d2+1)*func[1,d1,d2+1,d3]+(d3+1)*func[2,d1,d2,d3+1]
104 elif len(sh)==dim+2:
105 sh_out=(sh[0],)+sh[2:]
106 out=numarray.zeros(sh_out)
107 if dim==2:
108 for d1 in range(sh[2]-1):
109 for d2 in range(sh[3]-1):
110 for k in range(sh[0]):
111 out[k,d1,d2]=(d1+1)*func[k,0,d1+1,d2]+(d2+1)*func[k,1,d1,d2+1]
112 else:
113 for d1 in range(sh[2]-1):
114 for d2 in range(sh[3]-1):
115 for d3 in range(sh[4]-1):
116 for k in range(sh[0]):
117 out[k,d1,d2,d3]=(d1+1)*func[k,0,d1+1,d2,d3]+(d2+1)*func[k,1,d1,d2+1,d3]+(d3+1)*func[k,2,d1,d2,d3+1]
118 else:
119 raise RunTimeError,"sdfadsa"
120 return out
121
122 def multiplyByFandX(f,n,func,in_index,out_shape,out_index):
123 if n==0:
124 n1=1
125 else:
126 n1=0
127 if n==1:
128 n2=1
129 else:
130 n2=0
131 if n==2:
132 n3=1
133 else:
134 n3=0
135
136 sh=func.shape
137 dim=len(sh)-len(in_index)
138 sh_out=out_shape+sh[len(in_index):]
139 out=numarray.zeros(sh_out)
140 if len(in_index)==1:
141 if len(out_index)==1:
142 if dim==2:
143 for d1 in range(n1,sh[1]):
144 for d2 in range(n2,sh[2]):
145 out[out_index[0],d1,d2]=f*func[in_index[0],d1-n1,d2-n2]
146
147 else:
148 for d1 in range(n1,sh[1]):
149 for d2 in range(n2,sh[2]):
150 for d3 in range(n3,sh[3]):
151 out[out_index[0],d1,d2,d3]=f*func[in_index[0],d1-n1,d2-n2,d3-n3]
152
153 elif len(out_index)==2:
154 if dim==2:
155 for d1 in range(n1,sh[1]):
156 for d2 in range(n2,sh[2]):
157 out[out_index[0],out_index[1],d1,d2]=f*func[in_index[0],d1-n1,d2-n2]
158 else:
159 for d1 in range(n1,sh[1]):
160 for d2 in range(n2,sh[2]):
161 for d3 in range(n3,sh[3]):
162 out[out_index[0],out_index[1],d1,d2,d3]=f*func[in_index[0],d1-n1,d2-n2,d3-n3]
163
164 else:
165 raise RunTimeError,"sdfadsa"
166
167 elif len(in_index)==2:
168 if len(out_index)==1:
169 if dim==2:
170 for d1 in range(n1,sh[2]):
171 for d2 in range(n2,sh[3]):
172 out[out_index[0],d1,d2]=f*func[in_index[0],in_index[1],d1-n1,d2-n2]
173 else:
174 for d1 in range(n1,sh[2]):
175 for d2 in range(n2,sh[3]):
176 for d3 in range(n3,sh[4]):
177 out[out_index[0],d1,d2,d3]=f*func[in_index[0],in_index[1],d1-n1,d2-n2,d3-n3]
178
179 elif len(out_index)==2:
180 if dim==2:
181 for d1 in range(n1,sh[2]):
182 for d2 in range(n2,sh[3]):
183 out[out_index[0],out_index[1],d1,d2]=f*func[in_index[0],in_index[1],d1-n1,d2-n2]
184 else:
185 for d1 in range(n1,sh[2]):
186 for d2 in range(n2,sh[3]):
187 for d3 in range(n3,sh[4]):
188 out[out_index[0],out_index[1],d1,d2,d3]=f*func[in_index[0],in_index[1],d1-n1,d2-n2,d3-n3]
189
190 else:
191 raise RunTimeError,"sdfadsa"
192 else:
193 raise RunTimeError,"sdfadsa"
194 return out
195
196 def powTextforX(c,e):
197 if e==0:
198 out=""
199 elif e==1:
200 out="x[%s]"%c
201 else:
202 out="x[%s]**%s"%(c,e)
203 return out
204 def multTextFactors(*txt):
205 out=""
206 for t in txt:
207 if len(t)>0:
208 if len(out)==0:
209 out=t
210 else:
211 out+="*"+t
212 return out
213
214 def addSumTermText(txt,f,fac):
215 if len(fac)>0:
216 fac2="*"+fac
217 if f>0:
218 if f==1:
219 if len(txt)==0:
220 out="%s"%fac
221 else:
222 out="+%s"%fac
223 else:
224 if len(txt)==0:
225 out="%s%s"%(f,fac2)
226 else:
227 out="+%s%s"%(f,fac2)
228 elif f<0:
229 if f==1:
230 out="-%s"%fac
231 else:
232 if len(txt)==0:
233 out="(-%s)%s"%(-f,fac2)
234 else:
235 out="-%s%s"%(-f,fac2)
236 else:
237 out=""
238 else:
239 if f>0:
240 if f==1:
241 if len(txt)==0:
242 out="1"
243 else:
244 out="+1"
245 else:
246 if len(txt)==0:
247 out="%s"%f
248 else:
249 out="+%s"%f
250 elif f<0:
251 if f==1:
252 out="-1"
253 else:
254 if len(txt)==0:
255 out="(-%s)"%(-f,)
256 else:
257 out="-%s"%(-f,)
258 else:
259 out=""
260
261 return txt+out
262
263 def makeTermText(func):
264 sh=func.shape
265 out=""
266 if len(sh)==2:
267 for d1 in range(sh[0]):
268 for d2 in range(sh[1]):
269 fac=multTextFactors(powTextforX(0,d1),powTextforX(1,d2))
270 out=addSumTermText(out,func[d1,d2],fac)
271 else:
272 for d1 in range(sh[0]):
273 for d2 in range(sh[1]):
274 for d3 in range(sh[2]):
275 fac=multTextFactors(powTextforX(0,d1),powTextforX(1,d2),powTextforX(2,d3))
276 out=addSumTermText(out,func[d1,d2,d3],fac)
277 if len(out)==0:
278 return "0"
279 else:
280 return out
281 def makeContactDerivativeText(func,dim,func_name):
282 sh=func.shape
283 out=""
284 if len(sh)==dim+2:
285 if sh[0]==1:
286 out2=""
287 for j in range(sh[1]):
288 t=makeTermText(func[0,j])
289 if not t == "0":
290 if len(out2)>0: out2+="+"
291 if j==0: out2+=t
292 if len(out2)==0:
293 out+=2*intend+"%s=%s\n"%(func_name,0)
294 else:
295 out+=2*intend+"%s=%s\n"%(func_name,out2)
296 else:
297 out+=2*intend+"%s=Data(0.,%s,FunctionOnContactZero(self.domain))\n"%(func_name,str(tuple(sh[:1])))
298 for i in range(sh[0]):
299 out2=""
300 for j in range(sh[1]):
301 t=makeTermText(func[i,j])
302 if not t == "0":
303 if len(out2)>0: out2+="+"
304 out2+="%s"%t
305 if len(out2)>0:
306 out+=2*intend+"%s[%s]=%s\n"%(func_name,i,out2)
307 else:
308 raise RunTimeError,"KKK"
309 return out
310
311 def makeNormalDerivativeText(func,dim,func_name,add_jump=False):
312 sh=func.shape
313 out=""
314 if add_jump:
315 extra=")*jump"
316 extra_pre="("
317 else:
318 extra=""
319 extra_pre=""
320
321 out=2*intend+"x_boundary=FunctionOnBoundary(self.domain).getX()\n"
322 if dim==3:
323 out+=2*intend+"n=%swhereZero(x_boundary[0] ,self.ABS_TOL)*numarray.array([-1., 0., 0.])"%extra_pre
324 out+="+whereZero(x_boundary[0]-1.,self.ABS_TOL)*numarray.array([ 1., 0., 0.])"
325 out+="+whereZero(x_boundary[1] ,self.ABS_TOL)*numarray.array([ 0.,-1., 0.])"
326 out+="+whereZero(x_boundary[1]-1.,self.ABS_TOL)*numarray.array([ 0., 1., 0.])"
327 out+="+whereZero(x_boundary[2] ,self.ABS_TOL)*numarray.array([ 0., 0.,-1.])"
328 out+="+whereZero(x_boundary[2]-1.,self.ABS_TOL)*numarray.array([ 0., 0., 1.])%s\n"%extra
329 else:
330 out+=2*intend+"n=%swhereZero(x_boundary[0] ,self.ABS_TOL)*numarray.array([-1., 0.])"%extra_pre
331 out+="+whereZero(x_boundary[0]-1.,self.ABS_TOL)*numarray.array([ 1., 0.])"
332 out+="+whereZero(x_boundary[1] ,self.ABS_TOL)*numarray.array([ 0.,-1.])"
333 out+="+whereZero(x_boundary[1]-1.,self.ABS_TOL)*numarray.array([ 0., 1.])%s\n"%extra
334 if len(sh)==dim+2:
335 if sh[0]==1:
336 out2=""
337 for j in range(sh[1]):
338 t=makeTermText(func[0,j])
339 if not t == "0":
340 if len(out2)>0: out2+="+"
341 out2+="n[%s]*(%s)"%(j,t)
342 if len(out2)==0:
343 out+=2*intend+"%s=%s\n"%(func_name,0)
344 else:
345 out+=2*intend+"%s=%s\n"%(func_name,out2)
346 else:
347 out+=2*intend+"%s=Data(0.,%s,FunctionOnBoundary(self.domain))\n"%(func_name,str(tuple(sh[:1])))
348 for i in range(sh[0]):
349 out2=""
350 for j in range(sh[1]):
351 t=makeTermText(func[i,j])
352 if not t == "0":
353 if len(out2)>0: out2+="+"
354 out2+="n[%s]*(%s)"%(j,t)
355 if len(out2)>0:
356 out+=2*intend+"%s[%s]=%s\n"%(func_name,i,out2)
357 else:
358 raise RunTimeError,"KKK"
359 return out
360
361 def makeFunctionText(func,dim,func_name,add_jump=False):
362 sh=func.shape
363 out=""
364 if add_jump:
365 extra=")*jump"
366 extra_pre="("
367 else:
368 extra=""
369 extra_pre=""
370
371 if len(sh)==dim:
372 out+=2*intend+"%s=%s%s%s\n"%(func_name,extra_pre,makeTermText(func),extra)
373 elif len(sh)==dim+1:
374 if sh[0]==1:
375 out+=2*intend+"%s=%s%s%s\n"%(func_name,extra_pre,makeTermText(func[0]),extra)
376 else:
377 out+=2*intend+"%s=Data(0.,%s,ContinuousFunction(self.domain))\n"%(func_name,str(tuple(sh[:1])))
378 for i in range(sh[0]):
379 t=makeTermText(func[i])
380 if not t=="0": out+=2*intend+"%s[%s]=%s%s%s\n"%(func_name,i,extra_pre,t,extra)
381 elif len(sh)==dim+2:
382 if sh[0]==1:
383 if sh[1]==1:
384 out+=2*intend+"%s=%s%s%s\n"%(func_name,extra_pre,makeTermText(func[0,0]),extra)
385 else:
386 out+=2*intend+"%s=Data(0.,%s,ContinuousFunction(self.domain))\n"%(func_name,str(tuple(sh[1:2])))
387 for j in range(sh[1]):
388 t=makeTermText(func[0,j])
389 if not t=="0": out+=2*intend+"%s[%s]=%s%s%s\n"%(func_name,j,extra_pre,t,extra)
390 else:
391 if sh[1]==1:
392 out+=2*intend+"%s=Data(0.,%s,ContinuousFunction(self.domain))\n"%(func_name,str(tuple(sh[:1])))
393 for i in range(sh[0]):
394 t=makeTermText(func[i,0])
395 if not t=="0": out+=2*intend+"%s[%s]=%s%s%s\n"%(func_name,i,extra_pre,t,extra)
396 else:
397 out+=2*intend+"%s=Data(0.,%s,ContinuousFunction(self.domain))\n"%(func_name,str(tuple(sh[:2])))
398 for i in range(sh[0]):
399 for j in range(sh[1]):
400 t=makeTermText(func[i,j])
401 if not t=="0": out+=2*intend+"%s[%s,%s]=%s%s%s\n"%(func_name,i,j,extra_pre,t,extra)
402 return out
403 for s in [1,2]:
404 for d in [2,3]:
405 for typ in ["Strong","Weak", "Contact"]:
406 for coffo in ["Full"]: # ["Full", "Reduced" ]:
407 if coffo=="Full":
408 func_i="Function"
409 else:
410 func_i="ReducedFunction"
411 for case in ["Const", "Vario" ]:
412 for solo in [1, 2 ]:
413 order=solo
414 if typ in ["Strong","Weak"]:
415 # coefficient A:
416 if s==1:
417 for i in range(d):
418 for j in range(d):
419 test_func=makeTestSolution(order,s,d)
420 grad_test_func=makeGradient(test_func,d)
421 body2=makeFunctionText(test_func,d,"u")
422 body2+=2*intend+"A_test=Data(0.,(%d,%d),%s(self.domain))\n"%(d,d,func_i)
423 if case == "Const" :
424 f=int(8*random.random())+1
425 body2+=2*intend+"A_test[%d,%d]=%s\n"%(i,j,f)
426 A_x_grad=multiplyByFandX(f,-1,grad_test_func,(0,j),(1,d),(0,i))
427 else:
428 body2+=2*intend+"A_test[%d,%d]=x[%s]\n"%(i,j,i)
429 A_x_grad=multiplyByFandX(1.,i,grad_test_func,(0,j),(1,d),(0,i))
430 if typ=="Weak":
431 body2+=makeFunctionText(A_x_grad,d,"X_test")
432 args="A=A_test, X=X_test"
433 else:
434 div_A_x_grad=makeDiv(A_x_grad,d)
435 body2+=makeFunctionText(-div_A_x_grad,d,"Y_test")
436 body2+=makeNormalDerivativeText(A_x_grad,d,"y_test")
437 args="A=A_test, Y=Y_test, y=y_test"
438
439 makeTitle(d,coffo,solo,s,"A",case,typ,body2,mark="%s%s"%(i,j),pdeargs=args)
440 else:
441 for p in range(s):
442 for i in range(d):
443 for q in range(s):
444 for j in range(d):
445 test_func=makeTestSolution(order,s,d)
446 grad_test_func=makeGradient(test_func,d)
447 body2=makeFunctionText(test_func,d,"u")
448 body2+=2*intend+"A_test=Data(0.,(%d,%d,%d,%d),%s(self.domain))\n"%(s,d,s,d,func_i)
449 if case == "Const" :
450 f=int(8*random.random())+1
451 body2+=2*intend+"A_test[%d,%d,%d,%d]=%s\n"%(p,i,q,j,f)
452 A_x_grad=multiplyByFandX(f,-1,grad_test_func,(q,j),(s,d),(p,i))
453 else:
454 body2+=2*intend+"A_test[%d,%d,%d,%d]=x[%s]\n"%(p,i,q,j,i)
455 A_x_grad=multiplyByFandX(1.,i,grad_test_func,(q,j),(s,d),(p,i))
456 if typ=="Weak":
457 body2+=makeFunctionText(A_x_grad,d,"X_test")
458 args="A=A_test, X=X_test"
459 else:
460 div_A_x_grad=makeDiv(A_x_grad,d)
461 body2+=makeFunctionText(-div_A_x_grad,d,"Y_test")
462 body2+=makeNormalDerivativeText(A_x_grad,d,"y_test")
463 args="A=A_test, Y=Y_test, y=y_test"
464 makeTitle(d,coffo,solo,s,"A",case,typ,body2,mark="%s%s%s%s"%(p,i,q,j),pdeargs=args)
465 # coefficient B:
466 if s==1:
467 for i in range(d):
468 test_func=makeTestSolution(order,s,d)
469 body2=makeFunctionText(test_func,d,"u",typ=="Contact")
470 body2+=2*intend+"B_test=Data(0.,(%d,),%s(self.domain))\n"%(d,func_i)
471 if case == "Const" :
472 f=int(8*random.random())+1
473 body2+=2*intend+"B_test[%d]=%s\n"%(i,f)
474 B=multiplyByFandX(f,-1,test_func,(0,),(1,d),(0,i))
475 else:
476 body2+=2*intend+"B_test[%d]=x[%i]\n"%(i,i)
477 B=multiplyByFandX(1.,i,test_func,(0,),(1,d),(0,i))
478 if typ=="Weak":
479 body2+=makeFunctionText(B,d,"X_test")
480 args="B=B_test, X=X_test"
481 elif typ=="Strong":
482 div_B=makeDiv(B,d)
483 body2+=makeFunctionText(-div_B,d,"Y_test")
484 body2+=makeNormalDerivativeText(B,d,"y_test")
485 args="B=B_test, Y=Y_test, y=y_test"
486 else:
487 div_B=makeDiv(B,d)
488 body2+=makeFunctionText(-div_B,d,"Y_test",typ=="Contact")
489 body2+=makeNormalDerivativeText(B,d,"y_test",typ=="Contact")
490 body2+=makeContactDerivativeText(-B,d,"y_contact_test")
491 args="B=B_test, Y=Y_test, y=y_test, y_contact=y_contact_test"
492 makeTitle(d,coffo,solo,s,"B",case,typ,body2,mark="%s"%(i),pdeargs=args,add_jump=(typ=="Contact"))
493 else:
494 for p in range(s):
495 for i in range(d):
496 for q in range(s):
497 test_func=makeTestSolution(order,s,d)
498 grad_test_func=makeGradient(test_func,d)
499 body2=makeFunctionText(test_func,d,"u",typ=="Contact")
500 body2+=2*intend+"B_test=Data(0.,(%d,%d,%d),%s(self.domain))\n"%(s,d,s,func_i)
501 if case == "Const" :
502 f=int(8*random.random())+1
503 body2+=2*intend+"B_test[%d,%d,%d]=%s\n"%(p,i,q,f)
504 B=multiplyByFandX(f,-1,test_func,(q,),(s,d),(p,i))
505 else:
506 body2+=2*intend+"B_test[%d,%d,%d]=x[%i]\n"%(p,i,q,i)
507 B=multiplyByFandX(1.,i,test_func,(q,),(s,d),(p,i))
508 if typ=="Weak":
509 body2+=makeFunctionText(B,d,"X_test")
510 args="B=B_test, X=X_test"
511 elif typ=="Strong":
512 div_B=makeDiv(B,d)
513 body2+=makeFunctionText(-div_B,d,"Y_test")
514 body2+=makeNormalDerivativeText(B,d,"y_test")
515 args="B=B_test, Y=Y_test, y=y_test"
516 else:
517 div_B=makeDiv(B,d)
518 body2+=makeFunctionText(-div_B,d,"Y_test",typ=="Contact")
519 body2+=makeNormalDerivativeText(B,d,"y_test",typ=="Contact")
520 body2+=makeContactDerivativeText(B,d,"y_contact_test")
521 args="B=B_test, Y=Y_test, y=y_test, y_contact=y_contact_test"
522 makeTitle(d,coffo,solo,s,"B",case,typ,body2,mark="%s%s%s"%(p,i,q),pdeargs=args,add_jump=(typ=="Contact"))
523 if typ=="Strong":
524 # coefficient C:
525 if s==1:
526 for j in range(d):
527 test_func=makeTestSolution(order,s,d)
528 grad_test_func=makeGradient(test_func,d)
529 body2=makeFunctionText(test_func,d,"u")
530 body2+=2*intend+"C_test=Data(0.,(%d,),%s(self.domain))\n"%(d,func_i)
531 if case == "Const" :
532 f=int(8*random.random())+1
533 body2+=2*intend+"C_test[%d]=%s\n"%(j,f)
534 C_x_grad=multiplyByFandX(f,-1,grad_test_func,(0,j),(1,),(0,))
535 else:
536 body2+=2*intend+"C_test[%d]=x[%i]\n"%(j,j)
537 C_x_grad=multiplyByFandX(1,j,grad_test_func,(0,j),(1,),(0,))
538 body2+=makeFunctionText(C_x_grad,d,"Y_test")
539 args="C=C_test, Y=Y_test"
540 makeTitle(d,coffo,solo,s,"C",case,typ,body2,mark="%s"%(j),pdeargs=args)
541 else:
542 for p in range(s):
543 for q in range(s):
544 for j in range(d):
545 test_func=makeTestSolution(order,s,d)
546 grad_test_func=makeGradient(test_func,d)
547 body2=makeFunctionText(test_func,d,"u")
548 body2+=2*intend+"C_test=Data(0.,(%d,%d,%d),%s(self.domain))\n"%(s,s,d,func_i)
549 if case == "Const" :
550 f=int(8*random.random())+1
551 body2+=2*intend+"C_test[%d,%d,%d]=%s\n"%(p,q,j,f)
552 C_x_grad=multiplyByFandX(f,-1,grad_test_func,(q,j),(s,),(p,))
553 else:
554 body2+=2*intend+"C_test[%d,%d,%d]=x[%i]\n"%(p,q,j,j)
555 C_x_grad=multiplyByFandX(1.,j,grad_test_func,(q,j),(s,),(p,))
556 args="C=C_test, Y=Y_test"
557 body2+=makeFunctionText(C_x_grad,d,"Y_test")
558 makeTitle(d,coffo,solo,s,"C",case,typ,body2,mark="%s%s%s"%(p,q,j),pdeargs=args)
559 # coefficient D:
560 if s==1:
561 if case == "Const" :
562 test_func=makeTestSolution(order,s,d)
563 else:
564 test_func=makeTestSolution(order-1,s,d)
565 body2=makeFunctionText(test_func,d,"u")
566 if case == "Const" :
567 f=int(8*random.random())+1
568 body2+=2*intend+"D_test=Data(%s,(),%s(self.domain))\n"%(f,func_i)
569 D=multiplyByFandX(f,-1,test_func,(0,),(1,),(0,))
570 else:
571 body2+=2*intend+"D_test=%s(self.domain).getX()[0]\n"%func_i
572 D=multiplyByFandX(1,0,test_func,(0,),(1,),(0,))
573 args="D=D_test, Y=Y_test"
574 body2+=makeFunctionText(D,d,"Y_test")
575 makeTitle(d,coffo,solo,s,"D",case,typ,body2,mark="",pdeargs=args)
576 else:
577 for p in range(s):
578 for q in range(s):
579 if case == "Const" :
580 test_func=makeTestSolution(order,s,d)
581 else:
582 test_func=makeTestSolution(order-1,s,d)
583 body2=makeFunctionText(test_func,d,"u")
584 body2+=2*intend+"D_test=Data(0.,(%d,%d),%s(self.domain))\n"%(s,s,func_i)
585 if case == "Const" :
586 f=int(8*random.random())+1
587 body2+=2*intend+"D_test[%d,%d]=%s\n"%(p,q,f)
588 D=multiplyByFandX(f,-1,test_func,(q,),(s,),(p,))
589 else:
590 body2+=2*intend+"D_test[%d,%d]=x[0]\n"%(p,q)
591 D=multiplyByFandX(1.,0,test_func,(q,),(s,),(p,))
592 args="D=D_test, Y=Y_test"
593 body2+=makeFunctionText(D,d,"Y_test")
594 makeTitle(d,coffo,solo,s,"D",case,typ,body2,mark="%s%s"%(p,q),pdeargs=args)
595 if coffo=="Full":
596 # coefficient d:
597 if typ == "Strong":
598 if s==1:
599 if case == "Const" :
600 test_func=makeTestSolution(order,s,d)
601 else:
602 test_func=makeTestSolution(order-1,s,d)
603 body2=makeFunctionText(test_func,d,"u")
604 if case == "Const" :
605 f=int(8*random.random())+1
606 body2+=2*intend+"d_test=Data(%s,(),FunctionOnBoundary(self.domain))\n"%(f,)
607 D=multiplyByFandX(f,-1,test_func,(0,),(1,),(0,))
608 else:
609 body2+=2*intend+"d_test=interpolate(x[0],FunctionOnBoundary(self.domain))\n"
610 D=multiplyByFandX(1.,0,test_func,(0,),(1,),(0,))
611 body2+=makeFunctionText(D,d,"y_test")
612 args="d=d_test, y=y_test"
613 makeTitle(d,coffo,solo,s,"d",case,typ,body2,mark="",pdeargs=args)
614 else:
615 for p in range(s):
616 for q in range(s):
617 if case == "Const" :
618 test_func=makeTestSolution(order,s,d)
619 else:
620 test_func=makeTestSolution(order-1,s,d)
621 body2=makeFunctionText(test_func,d,"u")
622 body2+=2*intend+"d_test=Data(0.,(%d,%d),FunctionOnBoundary(self.domain))\n"%(s,s)
623 if case == "Const" :
624 f=int(8*random.random())+1
625 body2+=2*intend+"d_test[%d,%d]=%s\n"%(p,q,f)
626 D=multiplyByFandX(f,-1,test_func,(q,),(s,),(p,))
627 else:
628 body2+=2*intend+"d_test[%d,%d]=x[0]\n"%(p,q)
629
630 D=multiplyByFandX(1.,0,test_func,(q,),(s,),(p,))
631 args="d=d_test, y=y_test"
632 body2+=makeFunctionText(D,d,"y_test")
633 makeTitle(d,coffo,solo,s,"d",case,typ,body2,mark="%s%s"%(p,q),pdeargs=args)
634 # coefficient d_contact:
635 if typ == "Contact":
636 if s==1:
637 if case == "Const" :
638 test_func=makeTestSolution(order,s,d)
639 else:
640 test_func=makeTestSolution(order-1,s,d)
641
642 body2=makeFunctionText(test_func,d,"u",add_jump=True)
643 if case == "Const" :
644 f=int(8*random.random())+1
645 body2+=2*intend+"d_contact_test=Data(%s,(),FunctionOnContactZero(self.domain))\n"%(f,)
646 D=multiplyByFandX(f,-1,test_func,(0,),(1,),(0,))
647 else:
648 body2+=2*intend+"d_contact_test=interpolate(x[0],FunctionOnContactZero(self.domain))\n"
649 D=multiplyByFandX(1.,0,test_func,(0,),(1,),(0,))
650 body2+=makeFunctionText(D,d,"y_contact_test")
651 args="d_contact=d_contact_test, y_contact=y_contact_test"
652 makeTitle(d,coffo,solo,s,"d_contact",case,typ,body2,mark="",pdeargs=args,add_jump=True)
653 else:
654 for p in range(s):
655 for q in range(s):
656 if case == "Const" :
657 test_func=makeTestSolution(order,s,d)
658 else:
659 test_func=makeTestSolution(order-1,s,d)
660
661 body2=makeFunctionText(test_func,d,"u",add_jump=True)
662 body2+=2*intend+"d_contact_test=Data(0.,(%d,%d),FunctionOnContactZero(self.domain))\n"%(s,s)
663 if case == "Const" :
664 f=int(8*random.random())+1
665 body2+=2*intend+"d_contact_test[%d,%d]=%s\n"%(p,q,f)
666 D=multiplyByFandX(f,-1,test_func,(q,),(s,),(p,))
667 else:
668 body2+=2*intend+"d_contact_test[%d,%d]=x[0]\n"%(p,q)
669 D=multiplyByFandX(1.,0,test_func,(q,),(s,),(p,))
670 body2+=makeFunctionText(D,d,"y_contact_test")
671 args="d_contact=d_contact_test, y_contact=y_contact_test"
672 makeTitle(d,coffo,solo,s,"d_contact",case,typ,body2,mark="%s%s"%(p,q),pdeargs=args,add_jump=True
673 )
674 #================================
675 print "import unittest"
676 print "import numarray"
677 print "from esys.escript import *"
678 print "from esys.finley import Rectangle,Brick"
679 #print "class Test_assemblage_2Do1_Strong(unittest.TestCase):"
680 #print t_prog["2Do1_strong"]
681 #print "class Test_assemblage_2Do2_Strong(unittest.TestCase):"
682 #print t_prog["2Do2_strong"]
683 #print "class Test_assemblage_3Do1_Strong(unittest.TestCase):"
684 #print t_prog["3Do1_strong"]
685 #print "class Test_assemblage_3Do2_Strong(unittest.TestCase):"
686 #print t_prog["3Do2_strong"]
687 #print "class Test_assemblage_2Do1_Weak(unittest.TestCase):"
688 #print t_prog["2Do1_weak"]
689 #print "class Test_assemblage_2Do2_Weak(unittest.TestCase):"
690 #print t_prog["2Do2_weak"]
691 #print "class Test_assemblage_3Do1_Weak(unittest.TestCase):"
692 #print t_prog["3Do1_weak"]
693 #print "class Test_assemblage_3Do2_Weak(unittest.TestCase):"
694 #print t_prog["3Do2_weak"]
695 print "class Test_assemblage_2Do1_Contact(unittest.TestCase):"
696 print t_prog["2Do1_contact"]
697 #print "class Test_assemblage_2Do2_Contact(unittest.TestCase):"
698 #print t_prog["2Do2_contact"]
699 #print "class Test_assemblage_3Do1_Contact(unittest.TestCase):"
700 #print t_prog["3Do1_contact"]
701 #print "class Test_assemblage_3Do2_contact(unittest.TestCase):"
702 #print t_prog["3Do2_contact"]
703
704
705 print "from esys.escript import *\n"
706 print "from esys.finley import *\n"
707 print "from esys.escript.linearPDEs import LinearPDE\n"
708 print "import numarray"
709 print "import unittest"
710
711 # print "class Test_Finley_assemblage_2Do1(Test_assemblage_2Do1_Strong,Test_assemblage_2Do1_Weak):"
712 # print " RES_TOL=1.e-7"
713 # print " ABS_TOL=1.e-8"
714 # print " def setUp(self):"
715 # print " self.domain =Rectangle(2,2,1)"
716 print "class Test_Finley_assemblage_2Do1_Contact(Test_assemblage_2Do1_Contact):"
717 print " RES_TOL=1.e-7"
718 print " ABS_TOL=1.e-8"
719 print " def setUp(self):"
720 print " d1 = Rectangle(n0=2,n1=1,l0=0.5,order=1,useElementsOnFace=True)"
721 print " x1 = ContinuousFunction(d1).getX()"
722 print " ContinuousFunction(d1).setTags(1,Scalar(1,ContinuousFunction(d1)))"
723 print " d2 = Rectangle(n0=2,n1=1,l0=0.5,order=1,useElementsOnFace=True)"
724 print " ContinuousFunction(d2).setTags(2,Scalar(1,ContinuousFunction(d2)))"
725 print " d2.setX(d2.getX()+[0.5,0.])"
726 print " self.domain = JoinFaces([d1,d2])"
727 print "suite = unittest.TestSuite()"
728 print "suite.addTest(unittest.makeSuite(Test_Finley_assemblage_2Do1_Contact))"
729 print "unittest.TextTestRunner(verbosity=2).run(suite)"

  ViewVC Help
Powered by ViewVC 1.1.26