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

Revision 147 - (show annotations)
Fri Aug 12 01:45:47 2005 UTC (15 years, 8 months ago) by jgs
Original Path: trunk/esys2/finley/test/python/PoissonSolverTest.py
File MIME type: text/x-python
File size: 3212 byte(s)
```erge of development branch dev-02 back to main trunk on 2005-08-12

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

## Properties

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