test/python/AssembleTest.py

# assemble test: 
#
# $Id$

__copyright__="""  Copyright (c) 2006 by ACcESS MNRF
                    http://www.access.edu.au
                Primary Business: Queensland, Australia"""
__license__="""Licensed under the Open Software License version 3.0
             http://www.opensource.org/licenses/osl-3.0.php"""
from esys.escript import *
from esys.escript.linearPDEs import *
from esys import finley
from math import *

global seed
num_elem=2 # number of elements in each spatial direction
num_equations=3 # number of equations
integration_order=-1 # order of the integration scheme
order_u=1 # solution order


seed=1
rand_factor=sqrt(2.)
#
#   The test solution is represented in the basis 
#
#   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
#                 or 
#   B=[1,x_1,x_2,x_1**2,x_2**2,....,x_1**order_u,x_2**order_u]   for dim=2
#
#   length of B is len_Basis.
#
#   any test solution u with numComp components is represented by the numComp x len_Basis matrix sol.
#   the actual solution is then given by u=matmul(sol,B). 
#

global maxError,coeff,total_maxError,total_coeff
total_maxError=0
total_coeff=""
maxError=0
coeff=""
D=()


def randomNum():
     global seed
     seed+=1
     s=rand_factor*seed
     return s-int(s)
    
def generateRandom(arg):
    if len(arg)==0:
       return randomNum()
    out=numarray.ones(arg,numarray.Float64) # *randomNum()
    return out

def algebraicGrad(u):
   out=numarray.zeros([u.shape[0],dim,len_Basis],numarray.Float64)
   for s in range(u.shape[0]):
     for i in range(dim):
        for k in range(len_Basis):
          h=0
          for j in range(len_Basis):
             h+=u[s,j]*D[i][k,j]
          out[s,i,k]=h
   return out

def algebraicDiv(u):
   out=numarray.zeros([u.shape[0],len_Basis],numarray.Float64)
   for s in range(u.shape[0]):
     for k in range(len_Basis):
        h=0
        for i in range(dim):
          for j in range(len_Basis):
             h+=u[s,i,j]*D[i][k,j]
        out[s,k]=h
   return out

def mult4(A,u):
      out=numarray.zeros([A.shape[0],A.shape[1],len_Basis],numarray.Float64)
      for s in range(A.shape[0]):
        for i in range(A.shape[1]):
           for k in range(len_Basis):
             h=0
             for t in range(A.shape[2]):
               for j in range(A.shape[3]):
                 h+=A[s,i,t,j]*u[t,j,k]
             out[s,i,k]=h
      return out

def mult3_2(A,u):
      out=numarray.zeros([A.shape[0],A.shape[1],len_Basis],numarray.Float64)
      for s in range(A.shape[0]):
        for i in range(A.shape[1]):
           for k in range(len_Basis):
             h=0
             for t in range(A.shape[2]):
                 h+=A[s,i,t]*u[t,k]
             out[s,i,k]=h
      return out

def mult3_1(A,u):
      out=numarray.zeros([A.shape[0],len_Basis],numarray.Float64)
      for s in range(A.shape[0]):
           for k in range(len_Basis):
             h=0
             for t in range(A.shape[1]):
               for j in range(A.shape[2]):
                 h+=A[s,t,j]*u[t,j,k]
             out[s,k]=h
      return out

def mult2(A,u):
      out=numarray.zeros([A.shape[0],len_Basis],numarray.Float64)
      for s in range(A.shape[0]):
         for k in range(len_Basis):
           h=0
           for t in range(A.shape[1]):
               h+=A[s,t]*u[t,k]
           out[s,k]=h
      return out


def evaluate(u,this):
    x=this.getX()
    if u.rank==2:
       out=Data(value=0,shape=(u.shape[0],),what=this,expand=True)
       for i0 in range(u.shape[0]):
          out[i0]=u[i0,0]
          for p in range(order_u):
             for j in range(dim):
               out[i0]+=u[i0,p*dim+j+1]*x[j]**(p+1)
    else:
       if u.shape[0]==1:
          out=Data(value=0,shape=(u.shape[1],),what=this,expand=True)
          for i1 in range(u.shape[1]):
             out[i1]=u[0,i1,0]
             for p in range(order_u):
                for j in range(dim):
                  out[i1]+=u[0,i1,p*dim+j+1]*x[j]**(p+1)

       elif u.shape[1]==1:
          out=Data(value=0,shape=(u.shape[0],),what=this,expand=True)
          for i0 in range(u.shape[0]):
             out[i0]=u[i0,0,0]
             for p in range(order_u):
                for j in range(dim):
                  out[i0]+=u[i0,0,p*dim+j+1]*x[j]**(p+1)
       else:
          out=Data(value=0,shape=(u.shape[0],u.shape[1]),what=this,expand=True)
          for i0 in range(u.shape[0]):
           for i1 in range(u.shape[1]):
             out[i0,i1]=u[i0,i1,0]
             for p in range(order_u):
                for j in range(dim):
                  out[i0,i1]+=u[i0,i1,p*dim+j+1]*x[j]**(p+1)
    return out

def checkSystem(text,operator,u,rhs):
    global maxError,coeff
    error=Lsup(operator*u-rhs)
    print "@@ "+text+" error: ",error
    if error>=maxError:
       maxError=error
       coeff=text
#
#
def TestSystem(numEqu,numComp,mydomain,reduce):
   elem=Function(mydomain)
   face_elem=FunctionOnBoundary(mydomain)
   nodes=ContinuousFunction(mydomain)
   nrml=face_elem.getNormal()
   #
   #   test solution:
   #
   u=generateRandom([numComp,len_Basis])
   # u=numarray.zeros([numComp,len_Basis])
   # u[0,0]=1
   # u[0,1]=1
   U=evaluate(u,nodes)
   gradu=algebraicGrad(u)
   #
   #   test A:
   #
   for p in range(numEqu):
      for q in range(numComp):
         for i in range(dim):
            for j in range(dim):

               # check div( A grad(u) ) = div( X )
               c_A=numarray.zeros([numEqu,dim,numComp,dim])
               c_A[p,i,q,j]=1
               if numEqu==1 and numComp==1:
                  c_A2=numarray.reshape(c_A,[dim,dim])
                  text="A[%d,%d]"%(i,j)
               else:
                  c_A2=c_A
                  text="A[%d,%d,%d,%d]"%(p,i,q,j)
               x=mult4(c_A,gradu)
               mypde1=LinearPDE(mydomain)
               mypde1.setValue(A=c_A2,X=evaluate(x,elem))
               mypde1.setReducedOrderForSolutionTo(reduce)
               checkSystem(text+" const with X",mypde1.getOperator(),U,mypde1.getRightHandSide())
               # check div( A grad(u) ) = Y
               y=-algebraicDiv(x)
               mypde2=LinearPDE(mydomain)
               mypde2.setValue(Y=evaluate(y,elem),y=matrixmult(evaluate(x,face_elem),nrml))
               mypde2.setReducedOrderForSolutionsTo(reduce)
               checkSystem(text+" const with Y",mypde1.getOperator(),U,mypde2.getRightHandSide())

            # check div( B u ) = div( X )
            c_B=numarray.zeros([numEqu,dim,numComp])
            c_B[p,i,q]=1
            if numEqu==1 and numComp==1:
               c_B2=numarray.reshape(c_B,[dim])
               text="B[%d]"%(i)
            else:
               c_B2=c_B
               text="B[%d,%d,%d]"%(p,i,q)
            x=mult3_2(c_B,u)
            mypde1=LinearPDE(mydomain)
            mypde1.setValue(B=c_B2,X=evaluate(x,elem))
            mypde1.setReducedOrderForSolutionsTo(reduce)
            checkSystem(text+" const with X",mypde1.getOperator(),U,mypde1.getRightHandSide())
            # check div( B u ) = Y
            y=-algebraicDiv(x)
            mypde2=LinearPDE(mydomain)
            mypde2.setValue(Y=evaluate(y,elem),y=matrixmult(evaluate(x,face_elem),nrml))
            mypde2.setReducedOrderForSolutionsTo(reduce)
            checkSystem(text+" const with Y",mypde1.getOperator(),U,mypde2.getRightHandSide())

            # check C grad(u) = Y
            c_C=numarray.zeros([numEqu,numComp,dim])
            c_C[p,q,i]=1
            if numEqu==1 and numComp==1:
               c_C2=numarray.reshape(c_C,[dim])
               text="C[%d]"%(i)
            else:
               c_C2=c_C
               text="C[%d,%d,%d]"%(p,q,i)
            y=mult3_1(c_C,gradu)
            mypde1=LinearPDE(mydomain)
            mypde1.setValue(C=c_C2,Y=evaluate(y,elem))
            mypde1.setReducedOrderForSolutionsTo(reduce)
            checkSystem(text+" const with Y",mypde1.getOperator(),U,mypde1.getRightHandSide())


         # check D u= Y
         c_D=numarray.zeros([numEqu,numComp])
         c_D[p,q]=1
         if numEqu==1 and numComp==1:
            c_D2=numarray.reshape(c_D,[1])
            text="D"
         else:
            c_D2=c_D
            text="D[%d,%d]"%(p,q)
         y=mult2(c_D,u)
         mypde1=LinearPDE(mydomain)
         mypde1.setValue(D=c_D2,Y=evaluate(y,elem))
         mypde1.setReducedOrderForSolutionsTo(reduce)
         checkSystem(text+" const with Y",mypde1.getOperator(),U,mypde1.getRightHandSide())

#  we start:
for dim in [2,3]:
    len_Basis=dim*order_u+1
    #
    #  the differential operators:
    #   
    #  D_j(matmul(sol,B))=matmul(sol,D_jB)=matmul(matmul(sol,D[j]),B)
    #
    D=()
    for i in range(dim):
      D=D+(numarray.zeros([len_Basis,len_Basis],numarray.Float64),)
      for j in range(order_u):
         if j==0:
           D[i][0,i+j+1]=1
         else:
           D[i][(j-1)*dim+i+1,j*dim+i+1]=j+1
    #
    #   generate mydomain:
    #
    for order in [1,2]:
     for onElements in [False,True]:
        if onElements==True:
            onElmtext=", with elements on faces"
        else:
            onElmtext=""
        if dim==3:
           mydomain=finley.Brick(num_elem,num_elem,num_elem,order=order,integrationOrder=integration_order,useElementsOnFace=onElements)
        else:
           mydomain=finley.Rectangle(num_elem,num_elem,order=order,integrationOrder=integration_order,useElementsOnFace=onElements)
        for reduce in [False,True]:
           if reduce==True:
              redtext=",reduced"
           else:
              redtext=""
           #
           #  and start the test process:
           #
           for numEqu in range(1,num_equations+1):
               print "@@@ Start testing assembling with dim=%d and %d equations (order=%d%s%s)"%(dim,numEqu,order,redtext,onElmtext)
               TestSystem(numEqu,numEqu,mydomain,reduce)
               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)
               if maxError>=total_maxError:
                  total_maxError=maxError
                  total_coeff=coeff+" with %d equations in %dD (order=%d%s%s)"%(numEqu,dim,order,redtext,onElmtext)
    
print "@@@@ end testing assemblage with maximal error %e for %s"%(total_maxError,total_coeff)


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_12,x_22,x_32,....,x_1order_u,x_2order_u,x_3order_u] for dim=3
28	# or
29	# B=[1,x_1,x_2,x_12,x_22,....,x_1order_u,x_2order_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,pdim+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,pdim+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,pdim+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,pdim+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,jdim+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
Name	Value
svn:eol-style	native
svn:keywords	Author Date Id Revision