/[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 102 - (hide annotations)
Wed Dec 15 07:08:39 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: 3152 byte(s)
*** empty log message ***

1 jgs 102 # $Id$
2    
3     import sys
4     import os
5     import unittest
6    
7     esys_root=os.getenv('ESYS_ROOT')
8     sys.path.append(esys_root+'/finley/lib')
9     sys.path.append(esys_root+'/escript/lib')
10     sys.path.append(esys_root+'/escript/py_src')
11    
12     from escript import *
13     from util import *
14     from LinearPDEs import *
15     import finley
16    
17     ne_list=[10,15,22,33,50,75]
18     height_list=[0.25,0.5,1.]
19    
20    
21     def getDomain(dim,ne,height):
22    
23     if dim==2:
24     ne1=int(ne*height+0.5)
25     mydomain=finley.Rectangle(n0=ne,n1=ne1,l1=height,order=2)
26     totne=ne1*ne
27     else:
28     ne2=int(ne*height+0.5)
29     mydomain=finley.Brick(n0=ne,n1=ne,n2=ne2,l2=height,order=2)
30     totne=ne2*ne*ne
31     print "%d -dimensional domain generated."%dim
32     print "height of the domain is ",height
33     print "total number of elements is ",totne
34     return mydomain
35    
36    
37     def Solve1(mydomain,height):
38     print "Fully constraint solution"
39     l=[1.,1.,1.]
40     l[mydomain.getDim()-1]=height
41     cf=ContinuousFunction(mydomain)
42     x=cf.getX()
43     #construct exact solution:
44     u_ex=Scalar(1.,cf)
45     for i in range(mydomain.getDim()):
46     u_ex*=x[i]*(x[i]-l[i])
47     #construct mask:
48     msk=Scalar(0.,cf)
49     for i in range(mydomain.getDim()):
50     msk+=x[i].whereZero()+(x[i]-l[i]).whereZero()
51     #construct right hand side
52     f=Scalar(0,cf)
53     for i in range(mydomain.getDim()):
54     f_p=Scalar(1,cf)
55     for j in range(mydomain.getDim()):
56     if i==j:
57     f_p*=-2.
58     else:
59     f_p*=x[j]*(x[j]-l[j])
60     f+=f_p
61    
62     mypde=Poisson(f=f,q=msk)
63     u=mypde.getSolution()
64     error=Lsup(u-u_ex)/Lsup(u_ex)
65     print "error = ",error
66     return error
67    
68     def Solve2(mydomain,height):
69     print "Partially constraint solution"
70     l=[1.,1.,1.]
71     l[mydomain.getDim()-1]=height
72     print l
73     cf=ContinuousFunction(mydomain)
74     x=cf.getX()
75     #construct exact solution:
76     u_ex=Scalar(1.,cf)
77     for i in range(mydomain.getDim()):
78     u_ex*=x[i]*(2*l[i]-x[i])
79     #construct mask:
80     msk=Scalar(0.,cf)
81     for i in range(mydomain.getDim()):
82     msk+=x[i].whereZero()
83     #construct right hand side
84     f=Scalar(0,cf)
85     for i in range(mydomain.getDim()):
86     f_p=Scalar(1,cf)
87     for j in range(mydomain.getDim()):
88     if i==j:
89     f_p*=2.
90     else:
91     f_p*=x[j]*(2*l[j]-x[j])
92     f+=f_p
93     mypde=Poisson(f=f,q=msk)
94     u=mypde.getSolution()
95     error=Lsup(u-u_ex)/Lsup(u_ex)
96     print "error = ",error
97     return error
98    
99    
100     error=0
101     for ne in ne_list:
102     for dim in [2,3]:
103     for height in height_list:
104     print "***************************************************************"
105     mydomain= getDomain(dim,ne,height)
106     print "---------------------------------------------------------------"
107     error=max(error,Solve1(mydomain,height))
108     print "---------------------------------------------------------------"
109     error=max(error,Solve2(mydomain,height))
110     print "***************************************************************"
111    
112     print "***************************************************************"
113     print "maximum error: ",error
114     print "***************************************************************"

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26