# Annotation of /trunk/esys2/finley/test/python/AssembleTest.py

Revision 82 - (hide annotations)
Tue Oct 26 06:53:54 2004 UTC (17 years, 8 months ago) by jgs
File MIME type: text/x-python
File size: 10383 byte(s)
```Initial revision

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

## Properties

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

 ViewVC Help Powered by ViewVC 1.1.26