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

Revision 102 - (hide annotations)
Wed Dec 15 07:08:39 2004 UTC (16 years, 4 months ago) by jgs
File MIME type: text/x-python
File size: 3152 byte(s)
*** empty log message ***

 1 jgs 102 # \$Id\$ 2 3 import sys 4 import os 5 import unittest 6 7 esys_root=os.getenv('ESYS_ROOT') 8 sys.path.append(esys_root+'/finley/lib') 9 sys.path.append(esys_root+'/escript/lib') 10 sys.path.append(esys_root+'/escript/py_src') 11 12 from escript import * 13 from util import * 14 from LinearPDEs import * 15 import finley 16 17 ne_list=[10,15,22,33,50,75] 18 height_list=[0.25,0.5,1.] 19 20 21 def getDomain(dim,ne,height): 22 23 if dim==2: 24 ne1=int(ne*height+0.5) 25 mydomain=finley.Rectangle(n0=ne,n1=ne1,l1=height,order=2) 26 totne=ne1*ne 27 else: 28 ne2=int(ne*height+0.5) 29 mydomain=finley.Brick(n0=ne,n1=ne,n2=ne2,l2=height,order=2) 30 totne=ne2*ne*ne 31 print "%d -dimensional domain generated."%dim 32 print "height of the domain is ",height 33 print "total number of elements is ",totne 34 return mydomain 35 36 37 def Solve1(mydomain,height): 38 print "Fully constraint solution" 39 l=[1.,1.,1.] 40 l[mydomain.getDim()-1]=height 41 cf=ContinuousFunction(mydomain) 42 x=cf.getX() 43 #construct exact solution: 44 u_ex=Scalar(1.,cf) 45 for i in range(mydomain.getDim()): 46 u_ex*=x[i]*(x[i]-l[i]) 47 #construct mask: 48 msk=Scalar(0.,cf) 49 for i in range(mydomain.getDim()): 50 msk+=x[i].whereZero()+(x[i]-l[i]).whereZero() 51 #construct right hand side 52 f=Scalar(0,cf) 53 for i in range(mydomain.getDim()): 54 f_p=Scalar(1,cf) 55 for j in range(mydomain.getDim()): 56 if i==j: 57 f_p*=-2. 58 else: 59 f_p*=x[j]*(x[j]-l[j]) 60 f+=f_p 61 62 mypde=Poisson(f=f,q=msk) 63 u=mypde.getSolution() 64 error=Lsup(u-u_ex)/Lsup(u_ex) 65 print "error = ",error 66 return error 67 68 def Solve2(mydomain,height): 69 print "Partially constraint solution" 70 l=[1.,1.,1.] 71 l[mydomain.getDim()-1]=height 72 print l 73 cf=ContinuousFunction(mydomain) 74 x=cf.getX() 75 #construct exact solution: 76 u_ex=Scalar(1.,cf) 77 for i in range(mydomain.getDim()): 78 u_ex*=x[i]*(2*l[i]-x[i]) 79 #construct mask: 80 msk=Scalar(0.,cf) 81 for i in range(mydomain.getDim()): 82 msk+=x[i].whereZero() 83 #construct right hand side 84 f=Scalar(0,cf) 85 for i in range(mydomain.getDim()): 86 f_p=Scalar(1,cf) 87 for j in range(mydomain.getDim()): 88 if i==j: 89 f_p*=2. 90 else: 91 f_p*=x[j]*(2*l[j]-x[j]) 92 f+=f_p 93 mypde=Poisson(f=f,q=msk) 94 u=mypde.getSolution() 95 error=Lsup(u-u_ex)/Lsup(u_ex) 96 print "error = ",error 97 return error 98 99 100 error=0 101 for ne in ne_list: 102 for dim in [2,3]: 103 for height in height_list: 104 print "***************************************************************" 105 mydomain= getDomain(dim,ne,height) 106 print "---------------------------------------------------------------" 107 error=max(error,Solve1(mydomain,height)) 108 print "---------------------------------------------------------------" 109 error=max(error,Solve2(mydomain,height)) 110 print "***************************************************************" 111 112 print "***************************************************************" 113 print "maximum error: ",error 114 print "***************************************************************"

## Properties

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