/[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 104 - (hide annotations)
Fri Dec 17 07:43:12 2004 UTC (16 years, 4 months ago) by jgs
Original Path: trunk/esys2/finley/test/python/PoissonSolverTest.py
File MIME type: text/x-python
File size: 2964 byte(s)
*** empty log message ***

1 jgs 102 # $Id$
2    
3 jgs 104 from esys.escript import *
4     from esys.linearPDEs import Poisson
5     import esys.finley as finley
6 jgs 102
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(f=f,q=msk)
53     u=mypde.getSolution()
54     error=Lsup(u-u_ex)/Lsup(u_ex)
55     print "error = ",error
56     return error
57    
58     def Solve2(mydomain,height):
59     print "Partially constraint solution"
60     l=[1.,1.,1.]
61     l[mydomain.getDim()-1]=height
62     print l
63     cf=ContinuousFunction(mydomain)
64     x=cf.getX()
65     #construct exact solution:
66     u_ex=Scalar(1.,cf)
67     for i in range(mydomain.getDim()):
68     u_ex*=x[i]*(2*l[i]-x[i])
69     #construct mask:
70     msk=Scalar(0.,cf)
71     for i in range(mydomain.getDim()):
72     msk+=x[i].whereZero()
73     #construct right hand side
74     f=Scalar(0,cf)
75     for i in range(mydomain.getDim()):
76     f_p=Scalar(1,cf)
77     for j in range(mydomain.getDim()):
78     if i==j:
79     f_p*=2.
80     else:
81     f_p*=x[j]*(2*l[j]-x[j])
82     f+=f_p
83     mypde=Poisson(f=f,q=msk)
84     u=mypde.getSolution()
85     error=Lsup(u-u_ex)/Lsup(u_ex)
86     print "error = ",error
87     return error
88    
89    
90     error=0
91     for ne in ne_list:
92     for dim in [2,3]:
93     for height in height_list:
94     print "***************************************************************"
95     mydomain= getDomain(dim,ne,height)
96     print "---------------------------------------------------------------"
97     error=max(error,Solve1(mydomain,height))
98     print "---------------------------------------------------------------"
99     error=max(error,Solve2(mydomain,height))
100     print "***************************************************************"
101    
102     print "***************************************************************"
103     print "maximum error: ",error
104     print "***************************************************************"

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26