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

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

  ViewVC Help
Powered by ViewVC 1.1.26