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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 2344 - (show annotations)
Mon Mar 30 02:13:58 2009 UTC (10 years, 4 months ago) by jfenwick
File MIME type: text/x-python
File size: 3638 byte(s)
Change __url__ to launchpad site

1
2 ########################################################
3 #
4 # Copyright (c) 2003-2008 by University of Queensland
5 # Earth Systems Science Computational Center (ESSCC)
6 # http://www.uq.edu.au/esscc
7 #
8 # Primary Business: Queensland, Australia
9 # Licensed under the Open Software License version 3.0
10 # http://www.opensource.org/licenses/osl-3.0.php
11 #
12 ########################################################
13
14 __copyright__="""Copyright (c) 2003-2008 by University of Queensland
15 Earth Systems Science Computational Center (ESSCC)
16 http://www.uq.edu.au/esscc
17 Primary Business: Queensland, Australia"""
18 __license__="""Licensed under the Open Software License version 3.0
19 http://www.opensource.org/licenses/osl-3.0.php"""
20 __url__="https://launchpad.net/escript-finley"
21
22 import sys
23 import os
24 import unittest
25
26 from esys.escript import *
27 from esys.escript.linearPDEs import *
28 from esys import finley
29
30 import numarray
31
32 Pi=numarray.pi
33 numElements=15
34 sml=0.1
35
36 #
37 # this test the assemblage of problems with periodic boundary conditions:
38 #
39
40 #
41 # test solution is u_ex=sin(2*Pi*n*x0)*...*sin(2*Pi*n*x_dim)
42 #
43 def TheTest(msh,constraints,reduce):
44 n=ContinuousFunction(msh)
45 x=n.getX()
46 # set a test problem:
47 u_ex=Scalar(1,what=n)
48 for i in range(msh.getDim()):
49 u_ex*=sin(2*Pi*x[i])
50 mypde=LinearPDE(msh)
51 mypde.setSymmetryOn()
52 mypde.setDebugOn()
53 mypde.setReducedOrderTo(reduce)
54 mypde.setValue(A=numarray.identity(msh.getDim()),D=sml,Y=(sml+4*Pi**2*msh.getDim())*u_ex,q=constraints,r=u_ex)
55 return Lsup(mypde.getSolution()-u_ex)/Lsup(u_ex)
56
57 max_error=0
58 max_text="none"
59 for onElements in [False,True]:
60 if onElements==True:
61 onElmtext=", with elements on faces"
62 else:
63 onElmtext=""
64 for order in [1,2]:
65 for reduce in [False,True]:
66 if reduce==True:
67 redtext=",reduced"
68 else:
69 redtext=""
70 for dim in [2,3]:
71 if dim==2:
72 for i0 in [True,True]:
73 for i1 in [True,True]:
74 msh=finley.Rectangle(numElements,numElements,order,periodic0=i0,periodic1=i1,useElementsOnFace=onElements)
75 n=ContinuousFunction(msh)
76 x=n.getX()
77 c=Scalar(0,what=n)
78 if i0==False:
79 c+=whereZero(x[0])+whereZero(x[0]-1.)
80 if i1==False:
81 c+=whereZero(x[1])+whereZero(x[1]-1.)
82 error=TheTest(msh,c,reduce)
83 text="Rectangle order = %d%s%s, periodic0= %d, periodic1= %d: error= %f"%(order,redtext,onElmtext,i0,i1,error)
84 print "@@ ",text
85 if error>max_error:
86 max_error=error
87 max_text=text
88 elif dim==3:
89 for i0 in [True,False]:
90 for i1 in [True,False]:
91 for i2 in [True,False]:
92 msh=finley.Brick(numElements,numElements,numElements,order,periodic0=i0,periodic1=i1,periodic2=i2,useElementsOnFace=onElements)
93 n=ContinuousFunction(msh)
94 x=n.getX()
95 c=Scalar(0,what=n)
96 if i0==False:
97 c+=whereZero(x[0])+whereZero(x[0]-1.)
98 if i1==False:
99 c+=whereZero(x[1])+whereZero(x[1]-1.)
100 if i2==False:
101 c+=whereZero(x[2])+whereZero(x[2]-1.)
102 error=TheTest(msh,c,reduce)
103 text="Brick order = %d%s%s, periodic0= %d, periodic1= %d, periodic2= %d: error= %f"%(order,redtext,onElmtext,i0,i1,i2,error)
104 print "@@ ",text
105 if error>max_error:
106 max_error=error
107 max_text=text
108
109 print "@@@@ maximum error for :",max_text

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26