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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 108 - (show annotations)
Thu Jan 27 06:21:59 2005 UTC (14 years, 10 months ago) by jgs
File MIME type: text/x-python
File size: 3022 byte(s)
*** empty log message ***

1 # $Id$
2
3 from esys.escript import *
4 from esys.linearPDEs import Poisson
5 import esys.finley as finley
6
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(mydomain)
53 mypde.setValue(f=f,q=msk)
54 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 mypde=Poisson(mydomain)
85 mypde.setValue(f=f,q=msk)
86 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