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

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

  ViewVC Help
Powered by ViewVC 1.1.26