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

Annotation of /trunk/finley/test/python/gentest

Parent Directory Parent Directory | Revision Log Revision Log


Revision 768 - (hide annotations)
Fri Jun 30 09:26:16 2006 UTC (13 years ago) by gross
File size: 31774 byte(s)
and the generation of tests for contact elements works now.
1 gross 589 #!/usr/bin/python
2 gross 766 import numarray
3     import random
4 gross 589 t_prog={}
5 gross 766 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 gross 589 t_prog["2Do1_contact"]=""
14     t_prog["2Do2_contact"]=""
15     t_prog["3Do1_contact"]=""
16     t_prog["3Do2_contact"]=""
17     intend=" "
18 gross 766 max_order=3
19 gross 589
20 gross 766 def makeTitle(d,coeffo,o,s,coef,case,typ,body="",mark="",pdeargs="",add_jump=False):
21 gross 589 if mark=="":
22 gross 766 name="test_assemblage_%sD_solO%s_coeffO%s_neq%s_%s_%s_type%s"%(d,o,coeffo,s,coef,case,typ)
23 gross 589 else:
24 gross 766 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 gross 589 else:
85 gross 766 raise RunTimeError,"sdfadsa"
86     return out
87 gross 589
88 gross 766 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 gross 768 out=2*intend+"n_contact=FunctionOnContactZero(self.domain).getNormal()\n"
284 gross 766 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 gross 768 if j==0:
291     if len(out2)>0: out2+="+"
292     out2+=t
293 gross 766 if len(out2)==0:
294     out+=2*intend+"%s=%s\n"%(func_name,0)
295     else:
296 gross 768 out+=2*intend+"%s=n_contact[0]*(%s)\n"%(func_name,out2)
297 gross 766 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 gross 768 if j==0:
305     if len(out2)>0: out2+="+"
306     out2+="%s"%t
307 gross 766 if len(out2)>0:
308 gross 768 out+=2*intend+"%s[%s]=n_contact[0]*(%s)\n"%(func_name,i,out2)
309 gross 766 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 gross 768 out+=2*intend+"%s=%s%s%s\n"%(func_name,extra_pre,out2,extra)
348 gross 766 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 gross 768 out+=2*intend+"%s[%s]=%s%s%s\n"%(func_name,i,extra_pre,out2,extra)
359 gross 766 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 gross 589 for d in [2,3]:
407 gross 766 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 gross 589 for i in range(d):
420     for j in range(d):
421 gross 766 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 gross 589 for p in range(s):
444     for i in range(d):
445     for q in range(s):
446     for j in range(d):
447 gross 766 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 gross 589 # coefficient B:
468     if s==1:
469     for i in range(d):
470 gross 766 test_func=makeTestSolution(order,s,d)
471 gross 768 body2=makeFunctionText(test_func,d,"u")
472 gross 766 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 gross 768 body2+=makeFunctionText(-div_B,d,"Y_test")
491     body2+=makeNormalDerivativeText(B,d,"y_test")
492 gross 766 body2+=makeContactDerivativeText(-B,d,"y_contact_test")
493     args="B=B_test, Y=Y_test, y=y_test, y_contact=y_contact_test"
494 gross 768 makeTitle(d,coffo,solo,s,"B",case,typ,body2,mark="%s"%(i),pdeargs=args)
495 gross 589 else:
496     for p in range(s):
497     for i in range(d):
498     for q in range(s):
499 gross 766 test_func=makeTestSolution(order,s,d)
500     grad_test_func=makeGradient(test_func,d)
501 gross 768 body2=makeFunctionText(test_func,d,"u")
502 gross 766 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 gross 768 body2+=makeFunctionText(-div_B,d,"Y_test")
521     body2+=makeNormalDerivativeText(B,d,"y_test")
522     body2+=makeContactDerivativeText(-B,d,"y_contact_test")
523 gross 766 args="B=B_test, Y=Y_test, y=y_test, y_contact=y_contact_test"
524 gross 768 makeTitle(d,coffo,solo,s,"B",case,typ,body2,mark="%s%s%s"%(p,i,q),pdeargs=args)
525 gross 766 if typ=="Strong":
526 gross 589 # coefficient C:
527     if s==1:
528     for j in range(d):
529 gross 766 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 gross 589 else:
544     for p in range(s):
545 gross 766 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 gross 589 # coefficient D:
562     if s==1:
563 gross 766 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 gross 589 else:
579     for p in range(s):
580     for q in range(s):
581 gross 766 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 gross 589 print "import unittest"
678     print "import numarray"
679     print "from esys.escript import *"
680     print "from esys.finley import Rectangle,Brick"
681 gross 766 #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 gross 589 print t_prog["2Do1_contact"]
699 gross 766 #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 gross 589
706 gross 766
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 gross 589 print " RES_TOL=1.e-7"
720 gross 766 print " ABS_TOL=1.e-8"
721 gross 589 print " def setUp(self):"
722 gross 766 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 gross 589 print "suite = unittest.TestSuite()"
730 gross 766 print "suite.addTest(unittest.makeSuite(Test_Finley_assemblage_2Do1_Contact))"
731 gross 589 print "unittest.TextTestRunner(verbosity=2).run(suite)"

  ViewVC Help
Powered by ViewVC 1.1.26