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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 108 - (hide annotations)
Thu Jan 27 06:21:59 2005 UTC (14 years, 11 months ago) by jgs
File MIME type: text/x-python
File size: 3022 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 jgs 108 mypde=Poisson(mydomain)
53     mypde.setValue(f=f,q=msk)
54 jgs 102 u=mypde.getSolution()
55     error=Lsup(u-u_ex)/Lsup(u_ex)
56     print "error = ",error
57     return error
58    
59     def Solve2(mydomain,height):
60     print "Partially constraint solution"
61     l=[1.,1.,1.]
62     l[mydomain.getDim()-1]=height
63     print l
64     cf=ContinuousFunction(mydomain)
65     x=cf.getX()
66     #construct exact solution:
67     u_ex=Scalar(1.,cf)
68     for i in range(mydomain.getDim()):
69     u_ex*=x[i]*(2*l[i]-x[i])
70     #construct mask:
71     msk=Scalar(0.,cf)
72     for i in range(mydomain.getDim()):
73     msk+=x[i].whereZero()
74     #construct right hand side
75     f=Scalar(0,cf)
76     for i in range(mydomain.getDim()):
77     f_p=Scalar(1,cf)
78     for j in range(mydomain.getDim()):
79     if i==j:
80     f_p*=2.
81     else:
82     f_p*=x[j]*(2*l[j]-x[j])
83     f+=f_p
84 jgs 108 mypde=Poisson(mydomain)
85     mypde.setValue(f=f,q=msk)
86 jgs 102 u=mypde.getSolution()
87     error=Lsup(u-u_ex)/Lsup(u_ex)
88     print "error = ",error
89     return error
90    
91    
92     error=0
93     for ne in ne_list:
94     for dim in [2,3]:
95     for height in height_list:
96     print "***************************************************************"
97     mydomain= getDomain(dim,ne,height)
98     print "---------------------------------------------------------------"
99     error=max(error,Solve1(mydomain,height))
100     print "---------------------------------------------------------------"
101     error=max(error,Solve2(mydomain,height))
102     print "***************************************************************"
103    
104     print "***************************************************************"
105     print "maximum error: ",error
106     print "***************************************************************"

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26