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

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

Parent Directory Parent Directory | Revision Log Revision Log

Revision 1072 - (hide annotations)
Thu Mar 29 06:44:30 2007 UTC (14 years ago) by gross
File size: 34755 byte(s) 
PDE assemblage for reduced integration order + tests added.
1 gross 589 #!/usr/bin/python
2 gross 766 import numarray
3     import random
4 gross 589 t_prog={}
5 gross 769 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 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 769 name="test_assemblage_%sD_solO%s_coeffO%s_NEqu%s_%s_%s_type%s"%(d,o,coeffo,s,coef,case,typ)
23 gross 589 else:
24 gross 769 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 gross 766 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 gross 773 if d1+d2<=order: out[k,d1,d2]=getRandom()
57 gross 766 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 gross 773 if d1+d2+d3<=order: out[k,d1,d2,d3]=getRandom()
64 gross 766 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 gross 589 else:
90 gross 766 raise RunTimeError,"sdfadsa"
91     return out
92 gross 589
93 gross 766 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 gross 768 out=2*intend+"n_contact=FunctionOnContactZero(self.domain).getNormal()\n"
289 gross 766 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 gross 768 if j==0:
296     if len(out2)>0: out2+="+"
297     out2+=t
298 gross 766 if len(out2)==0:
299     out+=2*intend+"%s=%s\n"%(func_name,0)
300     else:
301 gross 768 out+=2*intend+"%s=n_contact[0]*(%s)\n"%(func_name,out2)
302 gross 766 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 gross 768 if j==0:
310     if len(out2)>0: out2+="+"
311     out2+="%s"%t
312 gross 766 if len(out2)>0:
313 gross 768 out+=2*intend+"%s[%s]=n_contact[0]*(%s)\n"%(func_name,i,out2)
314 gross 766 else:
315     raise RunTimeError,"KKK"
316     return out
317    
318 gross 1072 def makeNormalDerivativeText(func,dim,func_name,func_pre,add_jump=False):
319 gross 766 sh=func.shape
320     out=""
321     if add_jump:
322     extra=")*jump"
323     extra_pre="("
324     else:
325     extra=""
326     extra_pre=""
327    
328 gross 1072 out=2*intend+"x_boundary=%sFunctionOnBoundary(self.domain).getX()\n"%func_pre
329 gross 766 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 gross 768 out+=2*intend+"%s=%s%s%s\n"%(func_name,extra_pre,out2,extra)
353 gross 766 else:
354 gross 1072 out+=2*intend+"%s=Data(0.,%s,%sFunctionOnBoundary(self.domain))\n"%(func_name,str(tuple(sh[:1])),func_pre)
355 gross 766 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 gross 768 out+=2*intend+"%s[%s]=%s%s%s\n"%(func_name,i,extra_pre,out2,extra)
364 gross 766 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 gross 769 for s in [1,2,3]:
411 gross 589 for d in [2,3]:
412 gross 766 for typ in ["Strong","Weak", "Contact"]:
413 gross 1072 # for coffo in ["Full", "Reduced" ]:
414     for coffo in [ "Reduced" ]:
415 gross 766 if coffo=="Full":
416 gross 1072 func_i=""
417     arg_post=""
418 gross 766 else:
419 gross 1072 func_i="Reduced"
420     arg_post="_reduced"
421 gross 766 for case in ["Const", "Vario" ]:
422     for solo in [1, 2 ]:
423     if typ in ["Strong","Weak"]:
424     # coefficient A:
425 gross 773 if case=="Vario" and typ=="Weak":
426     order=solo
427 gross 1072 elif coffo=="Reduced" and case=="Vario":
428     order=1
429 gross 773 else:
430     order=solo
431 gross 766 if s==1:
432 gross 589 for i in range(d):
433     for j in range(d):
434 gross 766 test_func=makeTestSolution(order,s,d)
435     grad_test_func=makeGradient(test_func,d)
436     body2=makeFunctionText(test_func,d,"u")
437 gross 1072 body2+=2*intend+"A_test=Data(0.,(%d,%d),%sFunction(self.domain))\n"%(d,d,func_i)
438 gross 766 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 gross 1072 args="A%s=A_test, X%s=X_test"%(arg_post,arg_post)
448 gross 766 else:
449     div_A_x_grad=makeDiv(A_x_grad,d)
450     body2+=makeFunctionText(-div_A_x_grad,d,"Y_test")
451 gross 1072 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 gross 766
458     makeTitle(d,coffo,solo,s,"A",case,typ,body2,mark="%s%s"%(i,j),pdeargs=args)
459     else:
460 gross 589 for p in range(s):
461     for i in range(d):
462     for q in range(s):
463     for j in range(d):
464 gross 766 test_func=makeTestSolution(order,s,d)
465     grad_test_func=makeGradient(test_func,d)
466     body2=makeFunctionText(test_func,d,"u")
467 gross 1072 body2+=2*intend+"A_test=Data(0.,(%d,%d,%d,%d),%sFunction(self.domain))\n"%(s,d,s,d,func_i)
468 gross 766 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 gross 1072 args="A%s=A_test, X%s=X_test"%(arg_post,arg_post)
478 gross 766 else:
479     div_A_x_grad=makeDiv(A_x_grad,d)
480     body2+=makeFunctionText(-div_A_x_grad,d,"Y_test")
481 gross 1072 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 gross 766 makeTitle(d,coffo,solo,s,"A",case,typ,body2,mark="%s%s%s%s"%(p,i,q,j),pdeargs=args)
488 gross 589 # coefficient B:
489 gross 769 if typ in ["Strong","Weak"] or coffo=="Full":
490 gross 773 if case=="Vario" and typ=="Weak":
491     order=solo-1
492 gross 1072 elif coffo=="Reduced" and case=="Vario":
493     order=0
494     elif coffo=="Reduced":
495     order=solo-1
496 gross 773 else:
497     order=solo
498 gross 769 if s==1:
499 gross 589 for i in range(d):
500 gross 769 test_func=makeTestSolution(order,s,d)
501     body2=makeFunctionText(test_func,d,"u")
502 gross 1072 body2+=2*intend+"B_test=Data(0.,(%d,),%sFunction(self.domain))\n"%(d,func_i)
503 gross 769 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 gross 1072 args="B%s=B_test, X%s=X_test"%(arg_post,arg_post)
513 gross 769 elif typ=="Strong":
514     div_B=makeDiv(B,d)
515     body2+=makeFunctionText(-div_B,d,"Y_test")
516 gross 1072 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 gross 769 else:
523     div_B=makeDiv(B,d)
524     body2+=makeFunctionText(-div_B,d,"Y_test")
525 gross 1072 body2+=makeNormalDerivativeText(B,d,"y_test",func_i)
526 gross 769 body2+=makeContactDerivativeText(-B,d,"y_contact_test")
527 gross 1072 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 gross 769 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 gross 1072 body2+=2*intend+"B_test=Data(0.,(%d,%d,%d),%sFunction(self.domain))\n"%(s,d,s,func_i)
537 gross 769 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 gross 1072 args="B%s=B_test, X%s=X_test"%(arg_post,arg_post)
547 gross 769 elif typ=="Strong":
548     div_B=makeDiv(B,d)
549     body2+=makeFunctionText(-div_B,d,"Y_test")
550 gross 1072 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 gross 769 else:
557     div_B=makeDiv(B,d)
558     body2+=makeFunctionText(-div_B,d,"Y_test")
559 gross 1072 body2+=makeNormalDerivativeText(B,d,"y_test",func_i)
560 gross 769 body2+=makeContactDerivativeText(-B,d,"y_contact_test")
561 gross 1072 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 gross 769 makeTitle(d,coffo,solo,s,"B",case,typ,body2,mark="%s%s%s"%(p,i,q),pdeargs=args)
563 gross 766 if typ=="Strong":
564 gross 589 # coefficient C:
565 gross 773 order=solo
566 gross 589 if s==1:
567     for j in range(d):
568 gross 766 test_func=makeTestSolution(order,s,d)
569     grad_test_func=makeGradient(test_func,d)
570     body2=makeFunctionText(test_func,d,"u")
571 gross 1072 body2+=2*intend+"C_test=Data(0.,(%d,),%sFunction(self.domain))\n"%(d,func_i)
572 gross 766 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 gross 1072 args="C%s=C_test, Y%s=Y_test"%(arg_post,arg_post)
581 gross 766 makeTitle(d,coffo,solo,s,"C",case,typ,body2,mark="%s"%(j),pdeargs=args)
582 gross 589 else:
583     for p in range(s):
584 gross 766 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 gross 1072 body2+=2*intend+"C_test=Data(0.,(%d,%d,%d),%sFunction(self.domain))\n"%(s,s,d,func_i)
590 gross 766 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 gross 1072 args="C%s=C_test, Y%s=Y_test"%(arg_post,arg_post)
598 gross 766 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 gross 589 # coefficient D:
601 gross 773 if case=="Vario":
602     order=solo-1
603     else:
604     order=solo
605 gross 589 if s==1:
606 gross 773 test_func=makeTestSolution(order,s,d)
607 gross 766 body2=makeFunctionText(test_func,d,"u")
608     if case == "Const" :
609     f=int(8*random.random())+1
610 gross 1072 body2+=2*intend+"D_test=Data(%s,(),%sFunction(self.domain))\n"%(f,func_i)
611 gross 766 D=multiplyByFandX(f,-1,test_func,(0,),(1,),(0,))
612     else:
613 gross 1072 body2+=2*intend+"D_test=%sFunction(self.domain).getX()[0]\n"%func_i
614 gross 766 D=multiplyByFandX(1,0,test_func,(0,),(1,),(0,))
615 gross 1072 args="D%s=D_test, Y%s=Y_test"%(arg_post,arg_post)
616 gross 766 body2+=makeFunctionText(D,d,"Y_test")
617     makeTitle(d,coffo,solo,s,"D",case,typ,body2,mark="",pdeargs=args)
618 gross 589 else:
619     for p in range(s):
620     for q in range(s):
621 gross 773 test_func=makeTestSolution(order,s,d)
622 gross 766 body2=makeFunctionText(test_func,d,"u")
623 gross 1072 body2+=2*intend+"D_test=Data(0.,(%d,%d),%sFunction(self.domain))\n"%(s,s,func_i)
624 gross 766 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 gross 1072 args="D%s=D_test, Y%s=Y_test"%(arg_post,arg_post)
632 gross 766 body2+=makeFunctionText(D,d,"Y_test")
633     makeTitle(d,coffo,solo,s,"D",case,typ,body2,mark="%s%s"%(p,q),pdeargs=args)
634 gross 1072 if coffo in ["Full", "Reduced"]:
635 gross 773 if case=="Vario":
636     order=solo-1
637     else:
638     order=solo
639 gross 766 # coefficient d:
640     if typ == "Strong":
641     if s==1:
642 gross 773 test_func=makeTestSolution(order,s,d)
643 gross 766 body2=makeFunctionText(test_func,d,"u")
644     if case == "Const" :
645     f=int(8*random.random())+1
646 gross 1072 body2+=2*intend+"d_test=Data(%s,(),%sFunctionOnBoundary(self.domain))\n"%(f,func_i)
647 gross 766 D=multiplyByFandX(f,-1,test_func,(0,),(1,),(0,))
648     else:
649 gross 1072 body2+=2*intend+"d_test=interpolate(x[0],%sFunctionOnBoundary(self.domain))\n"%(func_i,)
650 gross 766 D=multiplyByFandX(1.,0,test_func,(0,),(1,),(0,))
651     body2+=makeFunctionText(D,d,"y_test")
652 gross 1072 args="d%s=d_test, y%s=y_test"%(arg_post,arg_post)
653 gross 766 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 gross 773 test_func=makeTestSolution(order,s,d)
658 gross 766 body2=makeFunctionText(test_func,d,"u")
659 gross 1072 body2+=2*intend+"d_test=Data(0.,(%d,%d),%sFunctionOnBoundary(self.domain))\n"%(s,s,func_i)
660 gross 766 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 gross 1072 args="d%s=d_test, y%s=y_test"%(arg_post,arg_post)
669 gross 766 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 gross 773 if case=="Vario":
674     order=solo-1
675     else:
676     order=solo
677 gross 766 if s==1:
678 gross 773 test_func=makeTestSolution(order,s,d)
679 gross 766 body2=makeFunctionText(test_func,d,"u",add_jump=True)
680     if case == "Const" :
681     f=int(8*random.random())+1
682 gross 1072 body2+=2*intend+"d_contact_test=Data(%s,(),%sFunctionOnContactZero(self.domain))\n"%(f,func_i)
683 gross 766 D=multiplyByFandX(f,-1,test_func,(0,),(1,),(0,))
684     else:
685 gross 1072 body2+=2*intend+"d_contact_test=interpolate(x[0],%sFunctionOnContactZero(self.domain))\n"%func_i
686 gross 766 D=multiplyByFandX(1.,0,test_func,(0,),(1,),(0,))
687     body2+=makeFunctionText(D,d,"y_contact_test")
688 gross 1072 args="d_contact%s=d_contact_test, y_contact%s=y_contact_test"%(arg_post,arg_post)
689 gross 766 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 gross 773 test_func=makeTestSolution(order,s,d)
694 gross 766 body2=makeFunctionText(test_func,d,"u",add_jump=True)
695 gross 1072 body2+=2*intend+"d_contact_test=Data(0.,(%d,%d),%sFunctionOnContactZero(self.domain))\n"%(s,s,func_i)
696 gross 766 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 gross 1072 args="d_contact%s=d_contact_test, y_contact%s=y_contact_test"%(arg_post,arg_post)
705 gross 766 makeTitle(d,coffo,solo,s,"d_contact",case,typ,body2,mark="%s%s"%(p,q),pdeargs=args,add_jump=True
706     )
707     #================================
708 gross 589 print "import unittest"
709     print "import numarray"
710     print "from esys.escript import *"
711     print "from esys.finley import Rectangle,Brick"
712 gross 769 print "class Test_assemblage_2Do1(unittest.TestCase):"
713 gross 1072 if len(t_prog["2Do1"])>0:
714     print t_prog["2Do1"]
715     else:
716     print " pass"
717 gross 769 print "class Test_assemblage_2Do1_Reduced(unittest.TestCase):"
718 gross 1072 if len(t_prog["2Do1_reduced"])>0:
719     print t_prog["2Do1_reduced"]
720     else:
721     print " pass"
722 gross 766 print "class Test_assemblage_2Do1_Contact(unittest.TestCase):"
723 gross 1072 if len(t_prog["2Do1_contact"])>0:
724     print t_prog["2Do1_contact"]
725     else:
726     print " pass"
727 gross 769
728 gross 773 print "class Test_assemblage_2Do2(unittest.TestCase):"
729 gross 1072 if len(t_prog["2Do2"])>0:
730     print t_prog["2Do2"]
731     else:
732     print " pass"
733 gross 773 print "class Test_assemblage_2Do2_Reduced(unittest.TestCase):"
734 gross 1072 if len(t_prog["2Do2_reduced"])>0:
735     print t_prog["2Do2_reduced"]
736     else:
737     print " pass"
738 gross 773 print "class Test_assemblage_2Do2_Contact(unittest.TestCase):"
739 gross 1072 if len(t_prog["2Do2_contact"])>0:
740     print t_prog["2Do2_contact"]
741     else:
742     print " pass"
743 gross 769
744 gross 773 print "class Test_assemblage_3Do1(unittest.TestCase):"
745 gross 1072 if len(t_prog["3Do1"])>0:
746     print t_prog["3Do1"]
747     else:
748     print " pass"
749 gross 773 print "class Test_assemblage_3Do1_Reduced(unittest.TestCase):"
750 gross 1072 if len(t_prog["3Do1_reduced"])>0:
751     print t_prog["3Do1_reduced"]
752     else:
753     print " pass"
754 gross 773 print "class Test_assemblage_3Do1_Contact(unittest.TestCase):"
755 gross 1072 if len(t_prog["3Do1_contact"])>0:
756     print t_prog["3Do1_contact"]
757     else:
758     print " pass"
759 gross 769
760 gross 773 print "class Test_assemblage_3Do2(unittest.TestCase):"
761 gross 1072 if len(t_prog["3Do2"])>0:
762     print t_prog["3Do2"]
763     else:
764     print " pass"
765 gross 773 print "class Test_assemblage_3Do2_Reduced(unittest.TestCase):"
766 gross 1072 if len(t_prog["3Do2_reduced"])>0:
767     print t_prog["3Do2_reduced"]
768     else:
769     print " pass"
770 gross 773 print "class Test_assemblage_3Do2_contact(unittest.TestCase):"
771 gross 1072 if len(t_prog["3Do2_contact"])>0:
772     print t_prog["3Do2_contact"]
773     else:
774     print " pass"
775 gross 589
776 gross 769 print "from esys.escript import *"
777     print "from esys.finley import *"
778     print "from esys.escript.linearPDEs import LinearPDE"
779 gross 766 print "import numarray"
780     print "import unittest"
781 gross 769 print "NE=2 # number of element sin each spatial direction (must be even)"
782 gross 766
783 gross 1072 print "class Test_Finley_assemblage_2Do1(Test_assemblage_2Do2_Reduced):"
784 gross 769 print " RES_TOL=1.e-7"
785     print " ABS_TOL=1.e-8"
786     print " def setUp(self):"
787 gross 1072 print " self.domain =Rectangle(NE,NE,2)"
788     print "class Test_Finley_assemblage_2Do1_Contact(Test_assemblage_2Do2_Contact):"
789 gross 589 print " RES_TOL=1.e-7"
790 gross 766 print " ABS_TOL=1.e-8"
791 gross 589 print " def setUp(self):"
792 gross 1072 print " d1 = Rectangle(n0=int(NE/2),n1=NE,l0=0.5,order=2,useElementsOnFace=True)"
793 gross 766 print " x1 = ContinuousFunction(d1).getX()"
794     print " ContinuousFunction(d1).setTags(1,Scalar(1,ContinuousFunction(d1)))"
795 gross 1072 print " d2 = Rectangle(n0=int(NE/2),n1=NE,l0=0.5,order=2,useElementsOnFace=True)"
796 gross 766 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 gross 589 print "suite = unittest.TestSuite()"
800 gross 769 print "suite.addTest(unittest.makeSuite(Test_Finley_assemblage_2Do1))"
801 gross 766 print "suite.addTest(unittest.makeSuite(Test_Finley_assemblage_2Do1_Contact))"
802 gross 589 print "unittest.TextTestRunner(verbosity=2).run(suite)"
  ViewVC Help
Powered by ViewVC 1.1.26