/[escript]/trunk/finley/test/python/PoissonSolverTest.py
ViewVC logotype

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 147 - (hide 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 jgs 102 # $Id$
2    
3    
4 jgs 147 # 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 jgs 102 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 jgs 108 mypde=Poisson(mydomain)
57 jgs 147 mypde.setTolerance(1.e-10)
58 jgs 108 mypde.setValue(f=f,q=msk)
59 jgs 102 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 jgs 108 mypde=Poisson(mydomain)
90 jgs 147 mypde.setTolerance(1.e-10)
91 jgs 108 mypde.setValue(f=f,q=msk)
92 jgs 102 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 jgs 147 # for dim in [2]:
102 jgs 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

  ViewVC Help
Powered by ViewVC 1.1.26