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

Diff of /trunk/finley/test/python/gentest

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

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

Legend:
Removed from v.589  
changed lines
  Added in v.766

  ViewVC Help
Powered by ViewVC 1.1.26