/[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 773 by gross, Fri Jul 7 10:10:00 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"]=""
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,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=FunctionOnBoundary(self.domain).getX()\n"
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,FunctionOnBoundary(self.domain))\n"%(func_name,str(tuple(sh[:1])))
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"]:          if coffo=="Full":
415              # coefficient A:              func_i="Function"
416              if s==1:          else:
417                func_i="ReducedFunction"
418            for case in ["Const", "Vario" ]:
419              for solo in [1, 2 ]:
420                if typ in ["Strong","Weak"]:
421                   # coefficient A:
422                   if case=="Vario" and typ=="Weak":
423                       order=solo
424                   else:
425                       order=solo
426                   if s==1:
427                  for i in range(d):                  for i in range(d):
428                    for j in range(d):                    for j in range(d):
429                      makeTitle(d,o,s,"A",case,type,body="",mark="%s%s"%(i,j))                      test_func=makeTestSolution(order,s,d)
430              else:                      grad_test_func=makeGradient(test_func,d)
431                        body2=makeFunctionText(test_func,d,"u")
432                        body2+=2*intend+"A_test=Data(0.,(%d,%d),%s(self.domain))\n"%(d,d,func_i)
433                        if case == "Const" :
434                           f=int(8*random.random())+1
435                           body2+=2*intend+"A_test[%d,%d]=%s\n"%(i,j,f)
436                           A_x_grad=multiplyByFandX(f,-1,grad_test_func,(0,j),(1,d),(0,i))
437                        else:
438                           body2+=2*intend+"A_test[%d,%d]=x[%s]\n"%(i,j,i)  
439                           A_x_grad=multiplyByFandX(1.,i,grad_test_func,(0,j),(1,d),(0,i))
440                        if typ=="Weak":
441                            body2+=makeFunctionText(A_x_grad,d,"X_test")
442                            args="A=A_test, X=X_test"
443                        else:
444                            div_A_x_grad=makeDiv(A_x_grad,d)
445                            body2+=makeFunctionText(-div_A_x_grad,d,"Y_test")
446                            body2+=makeNormalDerivativeText(A_x_grad,d,"y_test")
447                            args="A=A_test, Y=Y_test, y=y_test"
448    
449                        makeTitle(d,coffo,solo,s,"A",case,typ,body2,mark="%s%s"%(i,j),pdeargs=args)
450                   else:
451                  for p in range(s):                  for p in range(s):
452                    for i in range(d):                    for i in range(d):
453                      for q in range(s):                      for q in range(s):
454                        for j in range(d):                        for j in range(d):
455                           makeTitle(d,o,s,"A",case,type,body="",mark="%s%s%s%s"%(p,i,q,j))                           test_func=makeTestSolution(order,s,d)
456                             grad_test_func=makeGradient(test_func,d)
457                             body2=makeFunctionText(test_func,d,"u")
458                             body2+=2*intend+"A_test=Data(0.,(%d,%d,%d,%d),%s(self.domain))\n"%(s,d,s,d,func_i)
459                             if case == "Const" :
460                                f=int(8*random.random())+1
461                                body2+=2*intend+"A_test[%d,%d,%d,%d]=%s\n"%(p,i,q,j,f)
462                                A_x_grad=multiplyByFandX(f,-1,grad_test_func,(q,j),(s,d),(p,i))                          
463                             else:
464                                body2+=2*intend+"A_test[%d,%d,%d,%d]=x[%s]\n"%(p,i,q,j,i)
465                                A_x_grad=multiplyByFandX(1.,i,grad_test_func,(q,j),(s,d),(p,i))
466                             if typ=="Weak":
467                                body2+=makeFunctionText(A_x_grad,d,"X_test")
468                                args="A=A_test, X=X_test"
469                             else:
470                                div_A_x_grad=makeDiv(A_x_grad,d)
471                                body2+=makeFunctionText(-div_A_x_grad,d,"Y_test")
472                                body2+=makeNormalDerivativeText(A_x_grad,d,"y_test")
473                                args="A=A_test, Y=Y_test, y=y_test"
474                             makeTitle(d,coffo,solo,s,"A",case,typ,body2,mark="%s%s%s%s"%(p,i,q,j),pdeargs=args)
475              # coefficient B:              # coefficient B:
476              if s==1:              if typ in ["Strong","Weak"] or coffo=="Full":
477                  for i in range(d):                if case=="Vario" and typ=="Weak":
478                    makeTitle(d,o,s,"B",case,type,body="",mark="%s"%(i))                     order=solo-1
479              else:                else:
480                  for p in range(s):                     order=solo
481                  if s==1:
482                    for i in range(d):                    for i in range(d):
483                      for q in range(s):                      test_func=makeTestSolution(order,s,d)
484                        makeTitle(d,o,s,"B",case,type,body="",mark="%s%s%s"%(p,i,q))                      body2=makeFunctionText(test_func,d,"u")
485              if type=="strong":                      body2+=2*intend+"B_test=Data(0.,(%d,),%s(self.domain))\n"%(d,func_i)
486                        if case == "Const" :
487                            f=int(8*random.random())+1
488                            body2+=2*intend+"B_test[%d]=%s\n"%(i,f)
489                            B=multiplyByFandX(f,-1,test_func,(0,),(1,d),(0,i))
490                        else:
491                            body2+=2*intend+"B_test[%d]=x[%i]\n"%(i,i)
492                            B=multiplyByFandX(1.,i,test_func,(0,),(1,d),(0,i))
493                        if typ=="Weak":
494                          body2+=makeFunctionText(B,d,"X_test")
495                          args="B=B_test, X=X_test"
496                        elif typ=="Strong":
497                          div_B=makeDiv(B,d)
498                          body2+=makeFunctionText(-div_B,d,"Y_test")
499                          body2+=makeNormalDerivativeText(B,d,"y_test")
500                          args="B=B_test, Y=Y_test, y=y_test"
501                        else:
502                          div_B=makeDiv(B,d)
503                          body2+=makeFunctionText(-div_B,d,"Y_test")
504                          body2+=makeNormalDerivativeText(B,d,"y_test")
505                          body2+=makeContactDerivativeText(-B,d,"y_contact_test")
506                          args="B=B_test, Y=Y_test, y=y_test, y_contact=y_contact_test"
507                        makeTitle(d,coffo,solo,s,"B",case,typ,body2,mark="%s"%(i),pdeargs=args)
508                  else:
509                      for p in range(s):
510                        for i in range(d):
511                          for q in range(s):
512                            test_func=makeTestSolution(order,s,d)      
513                            grad_test_func=makeGradient(test_func,d)
514                            body2=makeFunctionText(test_func,d,"u")
515                            body2+=2*intend+"B_test=Data(0.,(%d,%d,%d),%s(self.domain))\n"%(s,d,s,func_i)
516                            if case == "Const" :
517                               f=int(8*random.random())+1
518                               body2+=2*intend+"B_test[%d,%d,%d]=%s\n"%(p,i,q,f)
519                               B=multiplyByFandX(f,-1,test_func,(q,),(s,d),(p,i))
520                            else:
521                               body2+=2*intend+"B_test[%d,%d,%d]=x[%i]\n"%(p,i,q,i)
522                               B=multiplyByFandX(1.,i,test_func,(q,),(s,d),(p,i))
523                            if typ=="Weak":
524                               body2+=makeFunctionText(B,d,"X_test")
525                               args="B=B_test, X=X_test"
526                            elif typ=="Strong":
527                               div_B=makeDiv(B,d)
528                               body2+=makeFunctionText(-div_B,d,"Y_test")
529                               body2+=makeNormalDerivativeText(B,d,"y_test")
530                               args="B=B_test, Y=Y_test, y=y_test"
531                            else:
532                               div_B=makeDiv(B,d)
533                               body2+=makeFunctionText(-div_B,d,"Y_test")
534                               body2+=makeNormalDerivativeText(B,d,"y_test")
535                               body2+=makeContactDerivativeText(-B,d,"y_contact_test")
536                               args="B=B_test, Y=Y_test, y=y_test, y_contact=y_contact_test"
537                            makeTitle(d,coffo,solo,s,"B",case,typ,body2,mark="%s%s%s"%(p,i,q),pdeargs=args)
538                if typ=="Strong":
539                 # coefficient C:                 # coefficient C:
540                   order=solo
541                 if s==1:                 if s==1:
542                       for j in range(d):                       for j in range(d):
543                         makeTitle(d,o,s,"C",case,type,body="",mark="%s"%(j))                         test_func=makeTestSolution(order,s,d)
544                           grad_test_func=makeGradient(test_func,d)      
545                           body2=makeFunctionText(test_func,d,"u")
546                           body2+=2*intend+"C_test=Data(0.,(%d,),%s(self.domain))\n"%(d,func_i)
547                           if case == "Const" :
548                                 f=int(8*random.random())+1
549                                 body2+=2*intend+"C_test[%d]=%s\n"%(j,f)
550                                 C_x_grad=multiplyByFandX(f,-1,grad_test_func,(0,j),(1,),(0,))
551                           else:
552                                 body2+=2*intend+"C_test[%d]=x[%i]\n"%(j,j)
553                                 C_x_grad=multiplyByFandX(1,j,grad_test_func,(0,j),(1,),(0,))
554                           body2+=makeFunctionText(C_x_grad,d,"Y_test")
555                           args="C=C_test, Y=Y_test"
556                           makeTitle(d,coffo,solo,s,"C",case,typ,body2,mark="%s"%(j),pdeargs=args)
557                 else:                 else:
558                     for p in range(s):                     for p in range(s):
559                         for q in range(s):                       for q in range(s):
560                           for j in range(d):                         for j in range(d):
561                              makeTitle(d,o,s,"C",case,type,body="",mark="%s%s%s"%(p,q,j))                           test_func=makeTestSolution(order,s,d)      
562                             grad_test_func=makeGradient(test_func,d)      
563                             body2=makeFunctionText(test_func,d,"u")
564                             body2+=2*intend+"C_test=Data(0.,(%d,%d,%d),%s(self.domain))\n"%(s,s,d,func_i)
565                             if case == "Const" :
566                                 f=int(8*random.random())+1
567                                 body2+=2*intend+"C_test[%d,%d,%d]=%s\n"%(p,q,j,f)
568                                 C_x_grad=multiplyByFandX(f,-1,grad_test_func,(q,j),(s,),(p,))
569                             else:
570                                 body2+=2*intend+"C_test[%d,%d,%d]=x[%i]\n"%(p,q,j,j)
571                                 C_x_grad=multiplyByFandX(1.,j,grad_test_func,(q,j),(s,),(p,))
572                             args="C=C_test, Y=Y_test"
573                             body2+=makeFunctionText(C_x_grad,d,"Y_test")
574                             makeTitle(d,coffo,solo,s,"C",case,typ,body2,mark="%s%s%s"%(p,q,j),pdeargs=args)
575                 # coefficient D:                 # coefficient D:
576                   if case=="Vario":
577                       order=solo-1
578                   else:
579                       order=solo              
580                 if s==1:                 if s==1:
581                     makeTitle(d,o,s,"D",case,type,body="",mark="")                    test_func=makeTestSolution(order,s,d)      
582                      body2=makeFunctionText(test_func,d,"u")
583                      if case == "Const" :
584                          f=int(8*random.random())+1
585                          body2+=2*intend+"D_test=Data(%s,(),%s(self.domain))\n"%(f,func_i)
586                          D=multiplyByFandX(f,-1,test_func,(0,),(1,),(0,))
587                      else:
588                          body2+=2*intend+"D_test=%s(self.domain).getX()[0]\n"%func_i
589                          D=multiplyByFandX(1,0,test_func,(0,),(1,),(0,))
590                      args="D=D_test, Y=Y_test"
591                      body2+=makeFunctionText(D,d,"Y_test")
592                      makeTitle(d,coffo,solo,s,"D",case,typ,body2,mark="",pdeargs=args)
593                 else:                 else:
594                     for p in range(s):                     for p in range(s):
595                       for q in range(s):                       for q in range(s):
596                            makeTitle(d,o,s,"D",case,type,body="",mark="%s%s"%(p,q))                          test_func=makeTestSolution(order,s,d)      
597              # coefficient d:                          body2=makeFunctionText(test_func,d,"u")
598              if s==1:                          body2+=2*intend+"D_test=Data(0.,(%d,%d),%s(self.domain))\n"%(s,s,func_i)
599                   makeTitle(d,o,s,"d",case,type,body="",mark="")                          if case == "Const" :
600              else:                             f=int(8*random.random())+1
601                  for p in range(s):                             body2+=2*intend+"D_test[%d,%d]=%s\n"%(p,q,f)
602                    for q in range(s):                             D=multiplyByFandX(f,-1,test_func,(q,),(s,),(p,))
603                         makeTitle(d,o,s,"d_contact",case,type,body="",mark="%s%s"%(p,q))                          else:
604              # coefficient d_contact:                             body2+=2*intend+"D_test[%d,%d]=x[0]\n"%(p,q)
605              if s==1:                             D=multiplyByFandX(1.,0,test_func,(q,),(s,),(p,))
606                  makeTitle(d,o,s,"d_contact",case,type,body="",mark="")                          args="D=D_test, Y=Y_test"
607              else:                          body2+=makeFunctionText(D,d,"Y_test")                      
608                  for p in range(s):                          makeTitle(d,coffo,solo,s,"D",case,typ,body2,mark="%s%s"%(p,q),pdeargs=args)
609                    for q in range(s):              if coffo=="Full":
610                        makeTitle(d,o,s,"d_contact",case,type,body="",mark="%s%s"%(p,q))                 if case=="Vario":
611              #================================                     order=solo-1
612                   else:
613                       order=solo              
614                   # coefficient d:
615                   if typ == "Strong":
616                      if s==1:
617                         test_func=makeTestSolution(order,s,d)      
618                         body2=makeFunctionText(test_func,d,"u")
619                         if case == "Const" :
620                             f=int(8*random.random())+1
621                             body2+=2*intend+"d_test=Data(%s,(),FunctionOnBoundary(self.domain))\n"%(f,)
622                             D=multiplyByFandX(f,-1,test_func,(0,),(1,),(0,))
623                         else:
624                             body2+=2*intend+"d_test=interpolate(x[0],FunctionOnBoundary(self.domain))\n"
625                             D=multiplyByFandX(1.,0,test_func,(0,),(1,),(0,))
626                         body2+=makeFunctionText(D,d,"y_test")
627                         args="d=d_test, y=y_test"
628                         makeTitle(d,coffo,solo,s,"d",case,typ,body2,mark="",pdeargs=args)
629                      else:
630                        for p in range(s):
631                          for q in range(s):
632                            test_func=makeTestSolution(order,s,d)      
633                            body2=makeFunctionText(test_func,d,"u")
634                            body2+=2*intend+"d_test=Data(0.,(%d,%d),FunctionOnBoundary(self.domain))\n"%(s,s)
635                            if case == "Const" :
636                                  f=int(8*random.random())+1
637                                  body2+=2*intend+"d_test[%d,%d]=%s\n"%(p,q,f)
638                                  D=multiplyByFandX(f,-1,test_func,(q,),(s,),(p,))
639                            else:
640                                  body2+=2*intend+"d_test[%d,%d]=x[0]\n"%(p,q)
641                          
642                                  D=multiplyByFandX(1.,0,test_func,(q,),(s,),(p,))
643                            args="d=d_test, y=y_test"
644                            body2+=makeFunctionText(D,d,"y_test")
645                            makeTitle(d,coffo,solo,s,"d",case,typ,body2,mark="%s%s"%(p,q),pdeargs=args)
646                   # coefficient d_contact:
647                   if typ == "Contact":
648                      if case=="Vario":
649                         order=solo-1
650                      else:
651                         order=solo              
652                      if s==1:
653                        test_func=makeTestSolution(order,s,d)      
654                        body2=makeFunctionText(test_func,d,"u",add_jump=True)
655                        if case == "Const" :
656                             f=int(8*random.random())+1
657                             body2+=2*intend+"d_contact_test=Data(%s,(),FunctionOnContactZero(self.domain))\n"%(f,)
658                             D=multiplyByFandX(f,-1,test_func,(0,),(1,),(0,))
659                        else:
660                             body2+=2*intend+"d_contact_test=interpolate(x[0],FunctionOnContactZero(self.domain))\n"
661                             D=multiplyByFandX(1.,0,test_func,(0,),(1,),(0,))
662                        body2+=makeFunctionText(D,d,"y_contact_test")
663                        args="d_contact=d_contact_test, y_contact=y_contact_test"
664                        makeTitle(d,coffo,solo,s,"d_contact",case,typ,body2,mark="",pdeargs=args,add_jump=True)
665                      else:
666                       for p in range(s):
667                         for q in range(s):
668                            test_func=makeTestSolution(order,s,d)      
669                            body2=makeFunctionText(test_func,d,"u",add_jump=True)
670                            body2+=2*intend+"d_contact_test=Data(0.,(%d,%d),FunctionOnContactZero(self.domain))\n"%(s,s)
671                            if case == "Const" :
672                               f=int(8*random.random())+1
673                               body2+=2*intend+"d_contact_test[%d,%d]=%s\n"%(p,q,f)
674                               D=multiplyByFandX(f,-1,test_func,(q,),(s,),(p,))
675                            else:
676                               body2+=2*intend+"d_contact_test[%d,%d]=x[0]\n"%(p,q)
677                               D=multiplyByFandX(1.,0,test_func,(q,),(s,),(p,))
678                            body2+=makeFunctionText(D,d,"y_contact_test")
679                            args="d_contact=d_contact_test, y_contact=y_contact_test"
680                            makeTitle(d,coffo,solo,s,"d_contact",case,typ,body2,mark="%s%s"%(p,q),pdeargs=args,add_jump=True
681    )
682                   #================================
683  print "import unittest"  print "import unittest"
684  print "import numarray"  print "import numarray"
685  print "from esys.escript import *"  print "from esys.escript import *"
686  print "from esys.finley import Rectangle,Brick"  print "from esys.finley import Rectangle,Brick"
687  print "class Test_assemblage_2Do1(unittest.TestCase):"  print "class Test_assemblage_2Do1(unittest.TestCase):"
688  print t_prog["2Do1"]  print t_prog["2Do1"]
689    print "class Test_assemblage_2Do1_Reduced(unittest.TestCase):"
690    print t_prog["2Do1_reduced"]
691    print "class Test_assemblage_2Do1_Contact(unittest.TestCase):"
692    print t_prog["2Do1_contact"]
693    
694  print "class Test_assemblage_2Do2(unittest.TestCase):"  print "class Test_assemblage_2Do2(unittest.TestCase):"
695  print t_prog["2Do2"]  print t_prog["2Do2"]
696    print "class Test_assemblage_2Do2_Reduced(unittest.TestCase):"
697    print t_prog["2Do2_reduced"]
698    print "class Test_assemblage_2Do2_Contact(unittest.TestCase):"
699    print t_prog["2Do2_contact"]
700    
701  print "class Test_assemblage_3Do1(unittest.TestCase):"  print "class Test_assemblage_3Do1(unittest.TestCase):"
702  print t_prog["3Do1"]  print t_prog["3Do1"]
703    print "class Test_assemblage_3Do1_Reduced(unittest.TestCase):"
704    print t_prog["3Do1_reduced"]
705    print "class Test_assemblage_3Do1_Contact(unittest.TestCase):"
706    print t_prog["3Do1_contact"]
707    
708  print "class Test_assemblage_3Do2(unittest.TestCase):"  print "class Test_assemblage_3Do2(unittest.TestCase):"
709  print t_prog["3Do2"]  print t_prog["3Do2"]
710  print "class Test_assemblage_2Do1Contact(unittest.TestCase):"  print "class Test_assemblage_3Do2_Reduced(unittest.TestCase):"
711  print t_prog["2Do1_contact"]  print t_prog["3Do2_reduced"]
712  print "class Test_assemblage_2Do2Contact(unittest.TestCase):"  print "class Test_assemblage_3Do2_contact(unittest.TestCase):"
 print t_prog["2Do2_contact"]  
 print "class Test_assemblage_3Do1Contact(unittest.TestCase):"  
 print t_prog["3Do1_contact"]  
 print "class Test_assemblage_3Do2Contact(unittest.TestCase):"  
713  print t_prog["3Do2_contact"]  print t_prog["3Do2_contact"]
714    
715    
716    print "from esys.escript import *"
717    print "from esys.finley import *"
718    print "from esys.escript.linearPDEs import LinearPDE"
719    print "import numarray"
720    print "import unittest"
721    print "NE=2 # number of element sin each spatial direction (must be even)"
722    
723  print "class Test_Finley_assemblage_2Do1(Test_assemblage_2Do1):"  print "class Test_Finley_assemblage_2Do1(Test_assemblage_2Do1):"
724  print "   RES_TOL=1.e-7"  print "   RES_TOL=1.e-7"
725    print "   ABS_TOL=1.e-8"
726  print "   def setUp(self):"  print "   def setUp(self):"
727  print "       self.domain =Rectangle(2,2,1)"  print "       self.domain =Rectangle(NE,NE,1)"
728  print "class Test_Finley_assemblage_2Do2(Test_assemblage_2Do2):"  print "class Test_Finley_assemblage_2Do1_Contact(Test_assemblage_2Do1_Contact):"
 print "   RES_TOL=1.e-7"  
 print "   def setUp(self):"  
 print "       self.domain =Rectangle(2,2,2)"  
 print "class Test_Finley_assemblage_3Do1(Test_assemblage_3Do1):"  
 print "   RES_TOL=1.e-7"  
 print "   def setUp(self):"  
 print "       self.domain =Brick(2,2,2,1)"  
 print "class Test_Finley_assemblage_3Do2(Test_assemblage_3Do2):"  
729  print "   RES_TOL=1.e-7"  print "   RES_TOL=1.e-7"
730    print "   ABS_TOL=1.e-8"
731  print "   def setUp(self):"  print "   def setUp(self):"
732  print "       self.domain =Brick(2,2,2,2)"  print "       d1 = Rectangle(n0=int(NE/2),n1=NE,l0=0.5,order=1,useElementsOnFace=True)"
733    print "       x1 = ContinuousFunction(d1).getX()"
734    print "       ContinuousFunction(d1).setTags(1,Scalar(1,ContinuousFunction(d1)))"
735    print "       d2 = Rectangle(n0=int(NE/2),n1=NE,l0=0.5,order=1,useElementsOnFace=True)"
736    print "       ContinuousFunction(d2).setTags(2,Scalar(1,ContinuousFunction(d2)))"
737    print "       d2.setX(d2.getX()+[0.5,0.])"
738    print "       self.domain = JoinFaces([d1,d2])"
739  print "suite = unittest.TestSuite()"  print "suite = unittest.TestSuite()"
740  print "suite.addTest(unittest.makeSuite(Test_Finley_assemblage_2Do1))"  print "suite.addTest(unittest.makeSuite(Test_Finley_assemblage_2Do1))"
741  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))"  
742  print "unittest.TextTestRunner(verbosity=2).run(suite)"  print "unittest.TextTestRunner(verbosity=2).run(suite)"

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

  ViewVC Help
Powered by ViewVC 1.1.26