/[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 149 - (hide annotations)
Thu Sep 1 03:31:39 2005 UTC (15 years, 6 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

  ViewVC Help
Powered by ViewVC 1.1.26