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

Revision 149 - (show annotations)
Thu Sep 1 03:31:39 2005 UTC (14 years, 4 months ago) by jgs
File MIME type: text/x-python
File size: 10065 byte(s)
```Merge of development branch dev-02 back to main trunk on 2005-09-01

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

## Properties

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