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

Revision 108 - (show annotations)
Thu Jan 27 06:21:59 2005 UTC (15 years, 11 months ago) by jgs
File MIME type: text/x-python
File size: 3022 byte(s)
```*** empty log message ***

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

## Properties

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