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

Revision 149 - (hide annotations)
Thu Sep 1 03:31:39 2005 UTC (15 years, 7 months ago) by jgs
Original Path: trunk/esys2/finley/test/python/PoissonSolverTest.py
File MIME type: text/x-python
File size: 3109 byte(s)
```Merge of development branch dev-02 back to main trunk on 2005-09-01

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

## Properties

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