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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 768 - (show annotations)
Fri Jun 30 09:26:16 2006 UTC (12 years, 11 months ago) by gross
File size: 31774 byte(s)
and the generation of tests for contact elements works now.
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=2*intend+"n_contact=FunctionOnContactZero(self.domain).getNormal()\n"
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 j==0:
291 if len(out2)>0: out2+="+"
292 out2+=t
293 if len(out2)==0:
294 out+=2*intend+"%s=%s\n"%(func_name,0)
295 else:
296 out+=2*intend+"%s=n_contact[0]*(%s)\n"%(func_name,out2)
297 else:
298 out+=2*intend+"%s=Data(0.,%s,FunctionOnContactZero(self.domain))\n"%(func_name,str(tuple(sh[:1])))
299 for i in range(sh[0]):
300 out2=""
301 for j in range(sh[1]):
302 t=makeTermText(func[i,j])
303 if not t == "0":
304 if j==0:
305 if len(out2)>0: out2+="+"
306 out2+="%s"%t
307 if len(out2)>0:
308 out+=2*intend+"%s[%s]=n_contact[0]*(%s)\n"%(func_name,i,out2)
309 else:
310 raise RunTimeError,"KKK"
311 return out
312
313 def makeNormalDerivativeText(func,dim,func_name,add_jump=False):
314 sh=func.shape
315 out=""
316 if add_jump:
317 extra=")*jump"
318 extra_pre="("
319 else:
320 extra=""
321 extra_pre=""
322
323 out=2*intend+"x_boundary=FunctionOnBoundary(self.domain).getX()\n"
324 if dim==3:
325 out+=2*intend+"n=%swhereZero(x_boundary[0] ,self.ABS_TOL)*numarray.array([-1., 0., 0.])"%extra_pre
326 out+="+whereZero(x_boundary[0]-1.,self.ABS_TOL)*numarray.array([ 1., 0., 0.])"
327 out+="+whereZero(x_boundary[1] ,self.ABS_TOL)*numarray.array([ 0.,-1., 0.])"
328 out+="+whereZero(x_boundary[1]-1.,self.ABS_TOL)*numarray.array([ 0., 1., 0.])"
329 out+="+whereZero(x_boundary[2] ,self.ABS_TOL)*numarray.array([ 0., 0.,-1.])"
330 out+="+whereZero(x_boundary[2]-1.,self.ABS_TOL)*numarray.array([ 0., 0., 1.])%s\n"%extra
331 else:
332 out+=2*intend+"n=%swhereZero(x_boundary[0] ,self.ABS_TOL)*numarray.array([-1., 0.])"%extra_pre
333 out+="+whereZero(x_boundary[0]-1.,self.ABS_TOL)*numarray.array([ 1., 0.])"
334 out+="+whereZero(x_boundary[1] ,self.ABS_TOL)*numarray.array([ 0.,-1.])"
335 out+="+whereZero(x_boundary[1]-1.,self.ABS_TOL)*numarray.array([ 0., 1.])%s\n"%extra
336 if len(sh)==dim+2:
337 if sh[0]==1:
338 out2=""
339 for j in range(sh[1]):
340 t=makeTermText(func[0,j])
341 if not t == "0":
342 if len(out2)>0: out2+="+"
343 out2+="n[%s]*(%s)"%(j,t)
344 if len(out2)==0:
345 out+=2*intend+"%s=%s\n"%(func_name,0)
346 else:
347 out+=2*intend+"%s=%s%s%s\n"%(func_name,extra_pre,out2,extra)
348 else:
349 out+=2*intend+"%s=Data(0.,%s,FunctionOnBoundary(self.domain))\n"%(func_name,str(tuple(sh[:1])))
350 for i in range(sh[0]):
351 out2=""
352 for j in range(sh[1]):
353 t=makeTermText(func[i,j])
354 if not t == "0":
355 if len(out2)>0: out2+="+"
356 out2+="n[%s]*(%s)"%(j,t)
357 if len(out2)>0:
358 out+=2*intend+"%s[%s]=%s%s%s\n"%(func_name,i,extra_pre,out2,extra)
359 else:
360 raise RunTimeError,"KKK"
361 return out
362
363 def makeFunctionText(func,dim,func_name,add_jump=False):
364 sh=func.shape
365 out=""
366 if add_jump:
367 extra=")*jump"
368 extra_pre="("
369 else:
370 extra=""
371 extra_pre=""
372
373 if len(sh)==dim:
374 out+=2*intend+"%s=%s%s%s\n"%(func_name,extra_pre,makeTermText(func),extra)
375 elif len(sh)==dim+1:
376 if sh[0]==1:
377 out+=2*intend+"%s=%s%s%s\n"%(func_name,extra_pre,makeTermText(func[0]),extra)
378 else:
379 out+=2*intend+"%s=Data(0.,%s,ContinuousFunction(self.domain))\n"%(func_name,str(tuple(sh[:1])))
380 for i in range(sh[0]):
381 t=makeTermText(func[i])
382 if not t=="0": out+=2*intend+"%s[%s]=%s%s%s\n"%(func_name,i,extra_pre,t,extra)
383 elif len(sh)==dim+2:
384 if sh[0]==1:
385 if sh[1]==1:
386 out+=2*intend+"%s=%s%s%s\n"%(func_name,extra_pre,makeTermText(func[0,0]),extra)
387 else:
388 out+=2*intend+"%s=Data(0.,%s,ContinuousFunction(self.domain))\n"%(func_name,str(tuple(sh[1:2])))
389 for j in range(sh[1]):
390 t=makeTermText(func[0,j])
391 if not t=="0": out+=2*intend+"%s[%s]=%s%s%s\n"%(func_name,j,extra_pre,t,extra)
392 else:
393 if sh[1]==1:
394 out+=2*intend+"%s=Data(0.,%s,ContinuousFunction(self.domain))\n"%(func_name,str(tuple(sh[:1])))
395 for i in range(sh[0]):
396 t=makeTermText(func[i,0])
397 if not t=="0": out+=2*intend+"%s[%s]=%s%s%s\n"%(func_name,i,extra_pre,t,extra)
398 else:
399 out+=2*intend+"%s=Data(0.,%s,ContinuousFunction(self.domain))\n"%(func_name,str(tuple(sh[:2])))
400 for i in range(sh[0]):
401 for j in range(sh[1]):
402 t=makeTermText(func[i,j])
403 if not t=="0": out+=2*intend+"%s[%s,%s]=%s%s%s\n"%(func_name,i,j,extra_pre,t,extra)
404 return out
405 for s in [1,2]:
406 for d in [2,3]:
407 for typ in ["Strong","Weak", "Contact"]:
408 for coffo in ["Full"]: # ["Full", "Reduced" ]:
409 if coffo=="Full":
410 func_i="Function"
411 else:
412 func_i="ReducedFunction"
413 for case in ["Const", "Vario" ]:
414 for solo in [1, 2 ]:
415 order=solo
416 if typ in ["Strong","Weak"]:
417 # coefficient A:
418 if s==1:
419 for i in range(d):
420 for j in range(d):
421 test_func=makeTestSolution(order,s,d)
422 grad_test_func=makeGradient(test_func,d)
423 body2=makeFunctionText(test_func,d,"u")
424 body2+=2*intend+"A_test=Data(0.,(%d,%d),%s(self.domain))\n"%(d,d,func_i)
425 if case == "Const" :
426 f=int(8*random.random())+1
427 body2+=2*intend+"A_test[%d,%d]=%s\n"%(i,j,f)
428 A_x_grad=multiplyByFandX(f,-1,grad_test_func,(0,j),(1,d),(0,i))
429 else:
430 body2+=2*intend+"A_test[%d,%d]=x[%s]\n"%(i,j,i)
431 A_x_grad=multiplyByFandX(1.,i,grad_test_func,(0,j),(1,d),(0,i))
432 if typ=="Weak":
433 body2+=makeFunctionText(A_x_grad,d,"X_test")
434 args="A=A_test, X=X_test"
435 else:
436 div_A_x_grad=makeDiv(A_x_grad,d)
437 body2+=makeFunctionText(-div_A_x_grad,d,"Y_test")
438 body2+=makeNormalDerivativeText(A_x_grad,d,"y_test")
439 args="A=A_test, Y=Y_test, y=y_test"
440
441 makeTitle(d,coffo,solo,s,"A",case,typ,body2,mark="%s%s"%(i,j),pdeargs=args)
442 else:
443 for p in range(s):
444 for i in range(d):
445 for q in range(s):
446 for j in range(d):
447 test_func=makeTestSolution(order,s,d)
448 grad_test_func=makeGradient(test_func,d)
449 body2=makeFunctionText(test_func,d,"u")
450 body2+=2*intend+"A_test=Data(0.,(%d,%d,%d,%d),%s(self.domain))\n"%(s,d,s,d,func_i)
451 if case == "Const" :
452 f=int(8*random.random())+1
453 body2+=2*intend+"A_test[%d,%d,%d,%d]=%s\n"%(p,i,q,j,f)
454 A_x_grad=multiplyByFandX(f,-1,grad_test_func,(q,j),(s,d),(p,i))
455 else:
456 body2+=2*intend+"A_test[%d,%d,%d,%d]=x[%s]\n"%(p,i,q,j,i)
457 A_x_grad=multiplyByFandX(1.,i,grad_test_func,(q,j),(s,d),(p,i))
458 if typ=="Weak":
459 body2+=makeFunctionText(A_x_grad,d,"X_test")
460 args="A=A_test, X=X_test"
461 else:
462 div_A_x_grad=makeDiv(A_x_grad,d)
463 body2+=makeFunctionText(-div_A_x_grad,d,"Y_test")
464 body2+=makeNormalDerivativeText(A_x_grad,d,"y_test")
465 args="A=A_test, Y=Y_test, y=y_test"
466 makeTitle(d,coffo,solo,s,"A",case,typ,body2,mark="%s%s%s%s"%(p,i,q,j),pdeargs=args)
467 # coefficient B:
468 if s==1:
469 for i in range(d):
470 test_func=makeTestSolution(order,s,d)
471 body2=makeFunctionText(test_func,d,"u")
472 body2+=2*intend+"B_test=Data(0.,(%d,),%s(self.domain))\n"%(d,func_i)
473 if case == "Const" :
474 f=int(8*random.random())+1
475 body2+=2*intend+"B_test[%d]=%s\n"%(i,f)
476 B=multiplyByFandX(f,-1,test_func,(0,),(1,d),(0,i))
477 else:
478 body2+=2*intend+"B_test[%d]=x[%i]\n"%(i,i)
479 B=multiplyByFandX(1.,i,test_func,(0,),(1,d),(0,i))
480 if typ=="Weak":
481 body2+=makeFunctionText(B,d,"X_test")
482 args="B=B_test, X=X_test"
483 elif typ=="Strong":
484 div_B=makeDiv(B,d)
485 body2+=makeFunctionText(-div_B,d,"Y_test")
486 body2+=makeNormalDerivativeText(B,d,"y_test")
487 args="B=B_test, Y=Y_test, y=y_test"
488 else:
489 div_B=makeDiv(B,d)
490 body2+=makeFunctionText(-div_B,d,"Y_test")
491 body2+=makeNormalDerivativeText(B,d,"y_test")
492 body2+=makeContactDerivativeText(-B,d,"y_contact_test")
493 args="B=B_test, Y=Y_test, y=y_test, y_contact=y_contact_test"
494 makeTitle(d,coffo,solo,s,"B",case,typ,body2,mark="%s"%(i),pdeargs=args)
495 else:
496 for p in range(s):
497 for i in range(d):
498 for q in range(s):
499 test_func=makeTestSolution(order,s,d)
500 grad_test_func=makeGradient(test_func,d)
501 body2=makeFunctionText(test_func,d,"u")
502 body2+=2*intend+"B_test=Data(0.,(%d,%d,%d),%s(self.domain))\n"%(s,d,s,func_i)
503 if case == "Const" :
504 f=int(8*random.random())+1
505 body2+=2*intend+"B_test[%d,%d,%d]=%s\n"%(p,i,q,f)
506 B=multiplyByFandX(f,-1,test_func,(q,),(s,d),(p,i))
507 else:
508 body2+=2*intend+"B_test[%d,%d,%d]=x[%i]\n"%(p,i,q,i)
509 B=multiplyByFandX(1.,i,test_func,(q,),(s,d),(p,i))
510 if typ=="Weak":
511 body2+=makeFunctionText(B,d,"X_test")
512 args="B=B_test, X=X_test"
513 elif typ=="Strong":
514 div_B=makeDiv(B,d)
515 body2+=makeFunctionText(-div_B,d,"Y_test")
516 body2+=makeNormalDerivativeText(B,d,"y_test")
517 args="B=B_test, Y=Y_test, y=y_test"
518 else:
519 div_B=makeDiv(B,d)
520 body2+=makeFunctionText(-div_B,d,"Y_test")
521 body2+=makeNormalDerivativeText(B,d,"y_test")
522 body2+=makeContactDerivativeText(-B,d,"y_contact_test")
523 args="B=B_test, Y=Y_test, y=y_test, y_contact=y_contact_test"
524 makeTitle(d,coffo,solo,s,"B",case,typ,body2,mark="%s%s%s"%(p,i,q),pdeargs=args)
525 if typ=="Strong":
526 # coefficient C:
527 if s==1:
528 for j in range(d):
529 test_func=makeTestSolution(order,s,d)
530 grad_test_func=makeGradient(test_func,d)
531 body2=makeFunctionText(test_func,d,"u")
532 body2+=2*intend+"C_test=Data(0.,(%d,),%s(self.domain))\n"%(d,func_i)
533 if case == "Const" :
534 f=int(8*random.random())+1
535 body2+=2*intend+"C_test[%d]=%s\n"%(j,f)
536 C_x_grad=multiplyByFandX(f,-1,grad_test_func,(0,j),(1,),(0,))
537 else:
538 body2+=2*intend+"C_test[%d]=x[%i]\n"%(j,j)
539 C_x_grad=multiplyByFandX(1,j,grad_test_func,(0,j),(1,),(0,))
540 body2+=makeFunctionText(C_x_grad,d,"Y_test")
541 args="C=C_test, Y=Y_test"
542 makeTitle(d,coffo,solo,s,"C",case,typ,body2,mark="%s"%(j),pdeargs=args)
543 else:
544 for p in range(s):
545 for q in range(s):
546 for j in range(d):
547 test_func=makeTestSolution(order,s,d)
548 grad_test_func=makeGradient(test_func,d)
549 body2=makeFunctionText(test_func,d,"u")
550 body2+=2*intend+"C_test=Data(0.,(%d,%d,%d),%s(self.domain))\n"%(s,s,d,func_i)
551 if case == "Const" :
552 f=int(8*random.random())+1
553 body2+=2*intend+"C_test[%d,%d,%d]=%s\n"%(p,q,j,f)
554 C_x_grad=multiplyByFandX(f,-1,grad_test_func,(q,j),(s,),(p,))
555 else:
556 body2+=2*intend+"C_test[%d,%d,%d]=x[%i]\n"%(p,q,j,j)
557 C_x_grad=multiplyByFandX(1.,j,grad_test_func,(q,j),(s,),(p,))
558 args="C=C_test, Y=Y_test"
559 body2+=makeFunctionText(C_x_grad,d,"Y_test")
560 makeTitle(d,coffo,solo,s,"C",case,typ,body2,mark="%s%s%s"%(p,q,j),pdeargs=args)
561 # coefficient D:
562 if s==1:
563 if case == "Const" :
564 test_func=makeTestSolution(order,s,d)
565 else:
566 test_func=makeTestSolution(order-1,s,d)
567 body2=makeFunctionText(test_func,d,"u")
568 if case == "Const" :
569 f=int(8*random.random())+1
570 body2+=2*intend+"D_test=Data(%s,(),%s(self.domain))\n"%(f,func_i)
571 D=multiplyByFandX(f,-1,test_func,(0,),(1,),(0,))
572 else:
573 body2+=2*intend+"D_test=%s(self.domain).getX()[0]\n"%func_i
574 D=multiplyByFandX(1,0,test_func,(0,),(1,),(0,))
575 args="D=D_test, Y=Y_test"
576 body2+=makeFunctionText(D,d,"Y_test")
577 makeTitle(d,coffo,solo,s,"D",case,typ,body2,mark="",pdeargs=args)
578 else:
579 for p in range(s):
580 for q in range(s):
581 if case == "Const" :
582 test_func=makeTestSolution(order,s,d)
583 else:
584 test_func=makeTestSolution(order-1,s,d)
585 body2=makeFunctionText(test_func,d,"u")
586 body2+=2*intend+"D_test=Data(0.,(%d,%d),%s(self.domain))\n"%(s,s,func_i)
587 if case == "Const" :
588 f=int(8*random.random())+1
589 body2+=2*intend+"D_test[%d,%d]=%s\n"%(p,q,f)
590 D=multiplyByFandX(f,-1,test_func,(q,),(s,),(p,))
591 else:
592 body2+=2*intend+"D_test[%d,%d]=x[0]\n"%(p,q)
593 D=multiplyByFandX(1.,0,test_func,(q,),(s,),(p,))
594 args="D=D_test, Y=Y_test"
595 body2+=makeFunctionText(D,d,"Y_test")
596 makeTitle(d,coffo,solo,s,"D",case,typ,body2,mark="%s%s"%(p,q),pdeargs=args)
597 if coffo=="Full":
598 # coefficient d:
599 if typ == "Strong":
600 if s==1:
601 if case == "Const" :
602 test_func=makeTestSolution(order,s,d)
603 else:
604 test_func=makeTestSolution(order-1,s,d)
605 body2=makeFunctionText(test_func,d,"u")
606 if case == "Const" :
607 f=int(8*random.random())+1
608 body2+=2*intend+"d_test=Data(%s,(),FunctionOnBoundary(self.domain))\n"%(f,)
609 D=multiplyByFandX(f,-1,test_func,(0,),(1,),(0,))
610 else:
611 body2+=2*intend+"d_test=interpolate(x[0],FunctionOnBoundary(self.domain))\n"
612 D=multiplyByFandX(1.,0,test_func,(0,),(1,),(0,))
613 body2+=makeFunctionText(D,d,"y_test")
614 args="d=d_test, y=y_test"
615 makeTitle(d,coffo,solo,s,"d",case,typ,body2,mark="",pdeargs=args)
616 else:
617 for p in range(s):
618 for q in range(s):
619 if case == "Const" :
620 test_func=makeTestSolution(order,s,d)
621 else:
622 test_func=makeTestSolution(order-1,s,d)
623 body2=makeFunctionText(test_func,d,"u")
624 body2+=2*intend+"d_test=Data(0.,(%d,%d),FunctionOnBoundary(self.domain))\n"%(s,s)
625 if case == "Const" :
626 f=int(8*random.random())+1
627 body2+=2*intend+"d_test[%d,%d]=%s\n"%(p,q,f)
628 D=multiplyByFandX(f,-1,test_func,(q,),(s,),(p,))
629 else:
630 body2+=2*intend+"d_test[%d,%d]=x[0]\n"%(p,q)
631
632 D=multiplyByFandX(1.,0,test_func,(q,),(s,),(p,))
633 args="d=d_test, y=y_test"
634 body2+=makeFunctionText(D,d,"y_test")
635 makeTitle(d,coffo,solo,s,"d",case,typ,body2,mark="%s%s"%(p,q),pdeargs=args)
636 # coefficient d_contact:
637 if typ == "Contact":
638 if s==1:
639 if case == "Const" :
640 test_func=makeTestSolution(order,s,d)
641 else:
642 test_func=makeTestSolution(order-1,s,d)
643
644 body2=makeFunctionText(test_func,d,"u",add_jump=True)
645 if case == "Const" :
646 f=int(8*random.random())+1
647 body2+=2*intend+"d_contact_test=Data(%s,(),FunctionOnContactZero(self.domain))\n"%(f,)
648 D=multiplyByFandX(f,-1,test_func,(0,),(1,),(0,))
649 else:
650 body2+=2*intend+"d_contact_test=interpolate(x[0],FunctionOnContactZero(self.domain))\n"
651 D=multiplyByFandX(1.,0,test_func,(0,),(1,),(0,))
652 body2+=makeFunctionText(D,d,"y_contact_test")
653 args="d_contact=d_contact_test, y_contact=y_contact_test"
654 makeTitle(d,coffo,solo,s,"d_contact",case,typ,body2,mark="",pdeargs=args,add_jump=True)
655 else:
656 for p in range(s):
657 for q in range(s):
658 if case == "Const" :
659 test_func=makeTestSolution(order,s,d)
660 else:
661 test_func=makeTestSolution(order-1,s,d)
662
663 body2=makeFunctionText(test_func,d,"u",add_jump=True)
664 body2+=2*intend+"d_contact_test=Data(0.,(%d,%d),FunctionOnContactZero(self.domain))\n"%(s,s)
665 if case == "Const" :
666 f=int(8*random.random())+1
667 body2+=2*intend+"d_contact_test[%d,%d]=%s\n"%(p,q,f)
668 D=multiplyByFandX(f,-1,test_func,(q,),(s,),(p,))
669 else:
670 body2+=2*intend+"d_contact_test[%d,%d]=x[0]\n"%(p,q)
671 D=multiplyByFandX(1.,0,test_func,(q,),(s,),(p,))
672 body2+=makeFunctionText(D,d,"y_contact_test")
673 args="d_contact=d_contact_test, y_contact=y_contact_test"
674 makeTitle(d,coffo,solo,s,"d_contact",case,typ,body2,mark="%s%s"%(p,q),pdeargs=args,add_jump=True
675 )
676 #================================
677 print "import unittest"
678 print "import numarray"
679 print "from esys.escript import *"
680 print "from esys.finley import Rectangle,Brick"
681 #print "class Test_assemblage_2Do1_Strong(unittest.TestCase):"
682 #print t_prog["2Do1_strong"]
683 #print "class Test_assemblage_2Do2_Strong(unittest.TestCase):"
684 #print t_prog["2Do2_strong"]
685 #print "class Test_assemblage_3Do1_Strong(unittest.TestCase):"
686 #print t_prog["3Do1_strong"]
687 #print "class Test_assemblage_3Do2_Strong(unittest.TestCase):"
688 #print t_prog["3Do2_strong"]
689 #print "class Test_assemblage_2Do1_Weak(unittest.TestCase):"
690 #print t_prog["2Do1_weak"]
691 #print "class Test_assemblage_2Do2_Weak(unittest.TestCase):"
692 #print t_prog["2Do2_weak"]
693 #print "class Test_assemblage_3Do1_Weak(unittest.TestCase):"
694 #print t_prog["3Do1_weak"]
695 #print "class Test_assemblage_3Do2_Weak(unittest.TestCase):"
696 #print t_prog["3Do2_weak"]
697 print "class Test_assemblage_2Do1_Contact(unittest.TestCase):"
698 print t_prog["2Do1_contact"]
699 #print "class Test_assemblage_2Do2_Contact(unittest.TestCase):"
700 #print t_prog["2Do2_contact"]
701 #print "class Test_assemblage_3Do1_Contact(unittest.TestCase):"
702 #print t_prog["3Do1_contact"]
703 #print "class Test_assemblage_3Do2_contact(unittest.TestCase):"
704 #print t_prog["3Do2_contact"]
705
706
707 print "from esys.escript import *\n"
708 print "from esys.finley import *\n"
709 print "from esys.escript.linearPDEs import LinearPDE\n"
710 print "import numarray"
711 print "import unittest"
712
713 # print "class Test_Finley_assemblage_2Do1(Test_assemblage_2Do1_Strong,Test_assemblage_2Do1_Weak):"
714 # print " RES_TOL=1.e-7"
715 # print " ABS_TOL=1.e-8"
716 # print " def setUp(self):"
717 # print " self.domain =Rectangle(2,2,1)"
718 print "class Test_Finley_assemblage_2Do1_Contact(Test_assemblage_2Do1_Contact):"
719 print " RES_TOL=1.e-7"
720 print " ABS_TOL=1.e-8"
721 print " def setUp(self):"
722 print " d1 = Rectangle(n0=2,n1=1,l0=0.5,order=1,useElementsOnFace=True)"
723 print " x1 = ContinuousFunction(d1).getX()"
724 print " ContinuousFunction(d1).setTags(1,Scalar(1,ContinuousFunction(d1)))"
725 print " d2 = Rectangle(n0=2,n1=1,l0=0.5,order=1,useElementsOnFace=True)"
726 print " ContinuousFunction(d2).setTags(2,Scalar(1,ContinuousFunction(d2)))"
727 print " d2.setX(d2.getX()+[0.5,0.])"
728 print " self.domain = JoinFaces([d1,d2])"
729 print "suite = unittest.TestSuite()"
730 print "suite.addTest(unittest.makeSuite(Test_Finley_assemblage_2Do1_Contact))"
731 print "unittest.TextTestRunner(verbosity=2).run(suite)"

  ViewVC Help
Powered by ViewVC 1.1.26