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

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

Parent Directory Parent Directory | Revision Log Revision Log


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

  ViewVC Help
Powered by ViewVC 1.1.26