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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 147 - (show annotations)
Fri Aug 12 01:45:47 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: 3212 byte(s)
erge of development branch dev-02 back to main trunk on 2005-08-12

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