Revision 97 - (show annotations)
Tue Dec 14 05:39:33 2004 UTC (16 years, 11 months ago) by jgs
File MIME type: text/x-python
File size: 6668 byte(s)
```*** empty log message ***

```
 1 # \$Id\$ 2 """ 3 Test a grad, interpolate and integrate over the unit square. 4 5 The tests are very basic. 6 7 by Lutz Gross, ACcESS, University of Queensland, Australia, 2003. 8 Version \$Id\$ 9 """ 10 11 import sys 12 import os 13 import unittest 14 15 from esys.escript import * 16 from esys.linearPDEs import * 17 from esys import finley 18 19 from math import * 20 from numarray import array 21 22 numElements=10 23 max_error=0. 24 max_error_text="" 25 26 for dim in [2,3]: 27 for order in [1,2]: 28 for onFace in [True,False]: 29 30 print "\n" 31 print "-----------------------------------------------------------------------------------------------------" 32 print "dim: %d, order: %i, onFace: %s." % (dim, order, onFace) 33 print "-----------------------------------------------------------------------------------------------------" 34 35 if onFace: 36 onFaceText=", on elements" 37 else: 38 onFaceText="" 39 40 case="dim=%d, order=%d%s" % (dim,order,onFaceText) 41 42 if dim==2: 43 mydomain=finley.Rectangle(numElements,numElements,order=order,useElementsOnFace=onFace) 44 m00=[[1,0],[0,0]] 45 m01=[[0,1],[0,0]] 46 m11=[[0,0],[0,1]] 47 h=5 48 else: 49 mydomain=finley.Brick(numElements,numElements,numElements,order=order,useElementsOnFace=onFace) 50 m00=[[1,0,0],[0,0,0]] 51 m01=[[0,1,0],[0,0,0]] 52 m11=[[0,0,0],[0,1,0]] 53 h=7 54 55 n=ContinuousFunction(mydomain) 56 e=Function(mydomain) 57 f=FunctionOnBoundary(mydomain) 58 d=Solution(mydomain) 59 r=ReducedSolution(mydomain) 60 61 # 62 # test gradient 63 # 64 65 test="error gradient in interior (nodes)" 66 67 x=n.getX()[0:2] 68 g=grad(x**order+x[1]*[1,0]) 69 ref=order*x[0]**(order-1)*m00+m01+order*x[1]**(order-1)*m11 70 error_norm=Lsup(ref-g) 71 72 text="%s: %55s = %e" % (case, test, error_norm) 73 print "%s" % text 74 75 if error_norm>max_error: 76 max_error_text=text 77 max_error=error_norm 78 79 # 80 # test gradient on degrees of freedom 81 # 82 83 test="error gradient in interior (degrees of freedom)" 84 85 x=n.getX()[0:2].interpolate(d) 86 g=grad(x**order+x[1]*[1,0]) 87 ref=order*x[0]**(order-1)*m00+m01+order*x[1]**(order-1)*m11 88 error_norm=Lsup(ref-g)/Lsup(ref) 89 90 text="%s: %55s = %e" % (case, test, error_norm) 91 print "%s" % text 92 93 if error_norm>max_error: 94 max_error_text=text 95 max_error=error_norm 96 97 # 98 # test gradient on reduced degrees of freedom 99 # 100 101 test="error gradient in interior (reduced degrees of freedom)" 102 103 x=n.getX()[0:2].interpolate(r) 104 g=grad(x+x[1]*[1,0]) 105 ref=Scalar(1,what=r)*m00+m01+Scalar(1,what=r)*m11 106 error_norm=Lsup(ref-g)/Lsup(ref) 107 108 text="%s: %55s = %e" % (case, test, error_norm) 109 print "%s" % text 110 111 if error_norm>max_error: 112 max_error_text=text 113 max_error=error_norm 114 115 # 116 # test integration over volume 117 # 118 119 test="error volume integration" 120 121 x=e.getX()[0:2] 122 error=integrate(x**2+[0,2.]*x)-array([1./3.,1./3.+2*1./2.]) 123 error_norm=sqrt(numarray.innerproduct(error,error)) 124 125 text="%s: %55s = %e" % (case, test, error_norm) 126 print "%s" % text 127 128 if error_norm>max_error: 129 max_error_text=text 130 max_error=error_norm 131 132 if onFace: 133 134 # 135 # gradient on the boundary: 136 # 137 138 test="error gradient on boundary" 139 140 x=n.getX()[0:2] 141 g=grad(x**order+x[1]*[1,0],where=f) 142 x=f.getX()[0:2] 143 ref=order*x[0]**(order-1)*m00+m01+order*x[1]**(order-1)*m11 144 error_norm=Lsup(g-ref) 145 146 text="%s: %55s = %e" % (case, test, error_norm) 147 print "%s" % text 148 149 if error_norm>max_error: 150 max_error_text=text 151 max_error=error_norm 152 153 # 154 # test gradient on degrees of freedom 155 # 156 157 test="error gradient on boundary (degrees of freedom)" 158 159 x=n.getX()[0:2].interpolate(d) 160 g=grad(x**order+x[1]*[1,0],where=f) 161 x=f.getX()[0:2] 162 ref=order*x[0]**(order-1)*m00+m01+order*x[1]**(order-1)*m11 163 error_norm=Lsup(ref-g)/Lsup(ref) 164 165 text="%s: %55s = %e" % (case, test, error_norm) 166 print "%s" % text 167 168 if error_norm>max_error: 169 max_error_text=text 170 max_error=error_norm 171 172 # 173 # test gradient on reduced degrees of freedom 174 # 175 176 test="error gradient on boundary (reduced degrees of freedom)" 177 178 x=n.getX()[0:2].interpolate(r) 179 g=grad(x+x[1]*[1,0],where=f) 180 ref=Scalar(1,what=r)*m00+m01+Scalar(1,what=r)*m11 181 error_norm=Lsup(ref-g)/Lsup(ref) 182 183 text="%s: %55s = %e" % (case, test, error_norm) 184 print "%s" % text 185 186 if error_norm>max_error: 187 max_error_text=text 188 max_error=error_norm 189 190 # 191 # test integration over boundary 192 # 193 194 test="error boundary integration" 195 196 x=f.getX()[0:2] 197 error=integrate(x**2+[0,2.]*x)-array([h/3.,h/3.+2*(h-1)/2.]) 198 error_norm=sqrt(numarray.innerproduct(error,error)) 199 200 text="%s: %55s = %e" % (case, test, error_norm) 201 print "%s" % text 202 203 if error_norm>max_error: 204 max_error_text=text 205 max_error=error_norm 206 207 # 208 # normal test: 209 # 210 211 test="error normals" 212 213 refNormal=Vector(0,what=f) 214 if dim==3: 215 refNormal.setTaggedValue(2,[1,0,0]) 216 refNormal.setTaggedValue(1,[-1,0,0]) 217 refNormal.setTaggedValue(20,[0,1,0]) 218 refNormal.setTaggedValue(10,[0,-1,0]) 219 refNormal.setTaggedValue(200,[0,0,1]) 220 refNormal.setTaggedValue(100,[0,0,-1]) 221 else: 222 refNormal.setTaggedValue(2,[1,0]) 223 refNormal.setTaggedValue(1,[-1,0]) 224 refNormal.setTaggedValue(20,[0,1]) 225 refNormal.setTaggedValue(10,[0,-1]) 226 error_norm=Lsup(f.getNormal()-refNormal) 227 228 text="%s: %55s = %e" % (case, test, error_norm) 229 print "%s" % text 230 231 if error_norm>max_error: 232 max_error_text=text 233 max_error=error_norm 234 235 print "-----------------------------------------------------------------------------------------------------" 236 237 print "\n\n" 238 print "******************************************************************************************************************" 239 print "maximal error:", max_error_text 240 print "******************************************************************************************************************"

## Properties

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