/[escript]/trunk/finley/test/python/AssembleTest.py
ViewVC logotype

Contents of /trunk/finley/test/python/AssembleTest.py

Parent Directory Parent Directory | Revision Log Revision Log


Revision 617 - (show annotations)
Wed Mar 22 02:58:17 2006 UTC (13 years, 8 months ago) by elspeth
File MIME type: text/x-python
File size: 10353 byte(s)
More copyright.

1 # assemble test:
2 #
3 # $Id$
4
5 __copyright__=""" Copyright (c) 2006 by ACcESS MNRF
6 http://www.access.edu.au
7 Primary Business: Queensland, Australia"""
8 __license__="""Licensed under the Open Software License version 3.0
9 http://www.opensource.org/licenses/osl-3.0.php"""
10 from esys.escript import *
11 from esys.escript.linearPDEs import *
12 from esys import finley
13 from math import *
14
15 global seed
16 num_elem=2 # number of elements in each spatial direction
17 num_equations=3 # number of equations
18 integration_order=-1 # order of the integration scheme
19 order_u=1 # solution order
20
21
22 seed=1
23 rand_factor=sqrt(2.)
24 #
25 # The test solution is represented in the basis
26 #
27 # B=[1,x_1,x_2,x_3,x_1**2,x_2**2,x_3**2,....,x_1**order_u,x_2**order_u,x_3**order_u] for dim=3
28 # or
29 # B=[1,x_1,x_2,x_1**2,x_2**2,....,x_1**order_u,x_2**order_u] for dim=2
30 #
31 # length of B is len_Basis.
32 #
33 # any test solution u with numComp components is represented by the numComp x len_Basis matrix sol.
34 # the actual solution is then given by u=matmul(sol,B).
35 #
36
37 global maxError,coeff,total_maxError,total_coeff
38 total_maxError=0
39 total_coeff=""
40 maxError=0
41 coeff=""
42 D=()
43
44
45 def randomNum():
46 global seed
47 seed+=1
48 s=rand_factor*seed
49 return s-int(s)
50
51 def generateRandom(arg):
52 if len(arg)==0:
53 return randomNum()
54 out=numarray.ones(arg,numarray.Float64) # *randomNum()
55 return out
56
57 def algebraicGrad(u):
58 out=numarray.zeros([u.shape[0],dim,len_Basis],numarray.Float64)
59 for s in range(u.shape[0]):
60 for i in range(dim):
61 for k in range(len_Basis):
62 h=0
63 for j in range(len_Basis):
64 h+=u[s,j]*D[i][k,j]
65 out[s,i,k]=h
66 return out
67
68 def algebraicDiv(u):
69 out=numarray.zeros([u.shape[0],len_Basis],numarray.Float64)
70 for s in range(u.shape[0]):
71 for k in range(len_Basis):
72 h=0
73 for i in range(dim):
74 for j in range(len_Basis):
75 h+=u[s,i,j]*D[i][k,j]
76 out[s,k]=h
77 return out
78
79 def mult4(A,u):
80 out=numarray.zeros([A.shape[0],A.shape[1],len_Basis],numarray.Float64)
81 for s in range(A.shape[0]):
82 for i in range(A.shape[1]):
83 for k in range(len_Basis):
84 h=0
85 for t in range(A.shape[2]):
86 for j in range(A.shape[3]):
87 h+=A[s,i,t,j]*u[t,j,k]
88 out[s,i,k]=h
89 return out
90
91 def mult3_2(A,u):
92 out=numarray.zeros([A.shape[0],A.shape[1],len_Basis],numarray.Float64)
93 for s in range(A.shape[0]):
94 for i in range(A.shape[1]):
95 for k in range(len_Basis):
96 h=0
97 for t in range(A.shape[2]):
98 h+=A[s,i,t]*u[t,k]
99 out[s,i,k]=h
100 return out
101
102 def mult3_1(A,u):
103 out=numarray.zeros([A.shape[0],len_Basis],numarray.Float64)
104 for s in range(A.shape[0]):
105 for k in range(len_Basis):
106 h=0
107 for t in range(A.shape[1]):
108 for j in range(A.shape[2]):
109 h+=A[s,t,j]*u[t,j,k]
110 out[s,k]=h
111 return out
112
113 def mult2(A,u):
114 out=numarray.zeros([A.shape[0],len_Basis],numarray.Float64)
115 for s in range(A.shape[0]):
116 for k in range(len_Basis):
117 h=0
118 for t in range(A.shape[1]):
119 h+=A[s,t]*u[t,k]
120 out[s,k]=h
121 return out
122
123
124 def evaluate(u,this):
125 x=this.getX()
126 if u.rank==2:
127 out=Data(value=0,shape=(u.shape[0],),what=this,expand=True)
128 for i0 in range(u.shape[0]):
129 out[i0]=u[i0,0]
130 for p in range(order_u):
131 for j in range(dim):
132 out[i0]+=u[i0,p*dim+j+1]*x[j]**(p+1)
133 else:
134 if u.shape[0]==1:
135 out=Data(value=0,shape=(u.shape[1],),what=this,expand=True)
136 for i1 in range(u.shape[1]):
137 out[i1]=u[0,i1,0]
138 for p in range(order_u):
139 for j in range(dim):
140 out[i1]+=u[0,i1,p*dim+j+1]*x[j]**(p+1)
141
142 elif u.shape[1]==1:
143 out=Data(value=0,shape=(u.shape[0],),what=this,expand=True)
144 for i0 in range(u.shape[0]):
145 out[i0]=u[i0,0,0]
146 for p in range(order_u):
147 for j in range(dim):
148 out[i0]+=u[i0,0,p*dim+j+1]*x[j]**(p+1)
149 else:
150 out=Data(value=0,shape=(u.shape[0],u.shape[1]),what=this,expand=True)
151 for i0 in range(u.shape[0]):
152 for i1 in range(u.shape[1]):
153 out[i0,i1]=u[i0,i1,0]
154 for p in range(order_u):
155 for j in range(dim):
156 out[i0,i1]+=u[i0,i1,p*dim+j+1]*x[j]**(p+1)
157 return out
158
159 def checkSystem(text,operator,u,rhs):
160 global maxError,coeff
161 error=Lsup(operator*u-rhs)
162 print "@@ "+text+" error: ",error
163 if error>=maxError:
164 maxError=error
165 coeff=text
166 #
167 #
168 def TestSystem(numEqu,numComp,mydomain,reduce):
169 elem=Function(mydomain)
170 face_elem=FunctionOnBoundary(mydomain)
171 nodes=ContinuousFunction(mydomain)
172 nrml=face_elem.getNormal()
173 #
174 # test solution:
175 #
176 u=generateRandom([numComp,len_Basis])
177 # u=numarray.zeros([numComp,len_Basis])
178 # u[0,0]=1
179 # u[0,1]=1
180 U=evaluate(u,nodes)
181 gradu=algebraicGrad(u)
182 #
183 # test A:
184 #
185 for p in range(numEqu):
186 for q in range(numComp):
187 for i in range(dim):
188 for j in range(dim):
189
190 # check div( A grad(u) ) = div( X )
191 c_A=numarray.zeros([numEqu,dim,numComp,dim])
192 c_A[p,i,q,j]=1
193 if numEqu==1 and numComp==1:
194 c_A2=numarray.reshape(c_A,[dim,dim])
195 text="A[%d,%d]"%(i,j)
196 else:
197 c_A2=c_A
198 text="A[%d,%d,%d,%d]"%(p,i,q,j)
199 x=mult4(c_A,gradu)
200 mypde1=LinearPDE(mydomain)
201 mypde1.setValue(A=c_A2,X=evaluate(x,elem))
202 mypde1.setReducedOrderForSolutionTo(reduce)
203 checkSystem(text+" const with X",mypde1.getOperator(),U,mypde1.getRightHandSide())
204 # check div( A grad(u) ) = Y
205 y=-algebraicDiv(x)
206 mypde2=LinearPDE(mydomain)
207 mypde2.setValue(Y=evaluate(y,elem),y=matrixmult(evaluate(x,face_elem),nrml))
208 mypde2.setReducedOrderForSolutionsTo(reduce)
209 checkSystem(text+" const with Y",mypde1.getOperator(),U,mypde2.getRightHandSide())
210
211 # check div( B u ) = div( X )
212 c_B=numarray.zeros([numEqu,dim,numComp])
213 c_B[p,i,q]=1
214 if numEqu==1 and numComp==1:
215 c_B2=numarray.reshape(c_B,[dim])
216 text="B[%d]"%(i)
217 else:
218 c_B2=c_B
219 text="B[%d,%d,%d]"%(p,i,q)
220 x=mult3_2(c_B,u)
221 mypde1=LinearPDE(mydomain)
222 mypde1.setValue(B=c_B2,X=evaluate(x,elem))
223 mypde1.setReducedOrderForSolutionsTo(reduce)
224 checkSystem(text+" const with X",mypde1.getOperator(),U,mypde1.getRightHandSide())
225 # check div( B u ) = Y
226 y=-algebraicDiv(x)
227 mypde2=LinearPDE(mydomain)
228 mypde2.setValue(Y=evaluate(y,elem),y=matrixmult(evaluate(x,face_elem),nrml))
229 mypde2.setReducedOrderForSolutionsTo(reduce)
230 checkSystem(text+" const with Y",mypde1.getOperator(),U,mypde2.getRightHandSide())
231
232 # check C grad(u) = Y
233 c_C=numarray.zeros([numEqu,numComp,dim])
234 c_C[p,q,i]=1
235 if numEqu==1 and numComp==1:
236 c_C2=numarray.reshape(c_C,[dim])
237 text="C[%d]"%(i)
238 else:
239 c_C2=c_C
240 text="C[%d,%d,%d]"%(p,q,i)
241 y=mult3_1(c_C,gradu)
242 mypde1=LinearPDE(mydomain)
243 mypde1.setValue(C=c_C2,Y=evaluate(y,elem))
244 mypde1.setReducedOrderForSolutionsTo(reduce)
245 checkSystem(text+" const with Y",mypde1.getOperator(),U,mypde1.getRightHandSide())
246
247
248 # check D u= Y
249 c_D=numarray.zeros([numEqu,numComp])
250 c_D[p,q]=1
251 if numEqu==1 and numComp==1:
252 c_D2=numarray.reshape(c_D,[1])
253 text="D"
254 else:
255 c_D2=c_D
256 text="D[%d,%d]"%(p,q)
257 y=mult2(c_D,u)
258 mypde1=LinearPDE(mydomain)
259 mypde1.setValue(D=c_D2,Y=evaluate(y,elem))
260 mypde1.setReducedOrderForSolutionsTo(reduce)
261 checkSystem(text+" const with Y",mypde1.getOperator(),U,mypde1.getRightHandSide())
262
263 # we start:
264 for dim in [2,3]:
265 len_Basis=dim*order_u+1
266 #
267 # the differential operators:
268 #
269 # D_j(matmul(sol,B))=matmul(sol,D_jB)=matmul(matmul(sol,D[j]),B)
270 #
271 D=()
272 for i in range(dim):
273 D=D+(numarray.zeros([len_Basis,len_Basis],numarray.Float64),)
274 for j in range(order_u):
275 if j==0:
276 D[i][0,i+j+1]=1
277 else:
278 D[i][(j-1)*dim+i+1,j*dim+i+1]=j+1
279 #
280 # generate mydomain:
281 #
282 for order in [1,2]:
283 for onElements in [False,True]:
284 if onElements==True:
285 onElmtext=", with elements on faces"
286 else:
287 onElmtext=""
288 if dim==3:
289 mydomain=finley.Brick(num_elem,num_elem,num_elem,order=order,integrationOrder=integration_order,useElementsOnFace=onElements)
290 else:
291 mydomain=finley.Rectangle(num_elem,num_elem,order=order,integrationOrder=integration_order,useElementsOnFace=onElements)
292 for reduce in [False,True]:
293 if reduce==True:
294 redtext=",reduced"
295 else:
296 redtext=""
297 #
298 # and start the test process:
299 #
300 for numEqu in range(1,num_equations+1):
301 print "@@@ Start testing assembling with dim=%d and %d equations (order=%d%s%s)"%(dim,numEqu,order,redtext,onElmtext)
302 TestSystem(numEqu,numEqu,mydomain,reduce)
303 print "@@@ end testing assembling (order=%d%s%s) with %d equations in %dD with maximum error %e for %s"%(order,redtext,onElmtext,numEqu,dim,maxError,coeff)
304 if maxError>=total_maxError:
305 total_maxError=maxError
306 total_coeff=coeff+" with %d equations in %dD (order=%d%s%s)"%(numEqu,dim,order,redtext,onElmtext)
307
308 print "@@@@ end testing assemblage with maximal error %e for %s"%(total_maxError,total_coeff)
309
310

Properties

Name Value
svn:eol-style native
svn:keywords Author Date Id Revision

  ViewVC Help
Powered by ViewVC 1.1.26