1 |
jgs |
102 |
# $Id$ |
2 |
gross |
567 |
from esys.escript import * |
3 |
gross |
569 |
from esys.escript.linearPDEs import LinearPDE |
4 |
jgs |
108 |
from esys.finley import Rectangle |
5 |
jgs |
102 |
#... set some parameters ... |
6 |
jgs |
108 |
kappa=1. |
7 |
jgs |
102 |
omega=0.1 |
8 |
|
|
eta=10. |
9 |
|
|
#... generate domain ... |
10 |
jgs |
113 |
mydomain = Rectangle(l0=5.,l1=1.,n0=50, n1=10) |
11 |
jgs |
102 |
#... open PDE and set coefficients ... |
12 |
gross |
567 |
mypde=LinearPDE(mydomain) |
13 |
|
|
mypde.setSymmetryOn() |
14 |
jgs |
102 |
n=mydomain.getNormal() |
15 |
|
|
x=mydomain.getX() |
16 |
gross |
569 |
mypde.setValue(A=kappa*kronecker(mydomain),D=omega,Y=omega*x[0], \ |
17 |
|
|
d=eta,y=kappa*n[0]+eta*x[0]) |
18 |
jgs |
102 |
#... calculate error of the PDE solution ... |
19 |
|
|
u=mypde.getSolution() |
20 |
gross |
569 |
print "error is ",Lsup(u-x[0]) |
21 |
lkettle |
575 |
# output should be similar to "error is 1.e-7" |
22 |
|
|
saveVTK("x0.xml",sol=u) |
23 |
|
|
|