/[escript]/trunk/doc/examples/geotutorial/steadystate.py
ViewVC logotype

Diff of /trunk/doc/examples/geotutorial/steadystate.py

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

revision 5287 by jfenwick, Wed Apr 9 05:41:57 2014 UTC revision 5288 by sshaw, Tue Dec 2 23:18:40 2014 UTC
# Line 24  __url__="https://launchpad.net/escript-f Line 24  __url__="https://launchpad.net/escript-f
24  # import tools  # import tools
25  from esys.escript import *  from esys.escript import *
26  from esys.escript.linearPDEs import LinearPDE  from esys.escript.linearPDEs import LinearPDE
27  from esys.dudley import Rectangle  try:
28        from esys.dudley import Rectangle
29        HAVE_DUDLEY = True
30    except ImportError:
31        HAVE_DUDLEY = False
32  from esys.weipa import saveVTK  from esys.weipa import saveVTK
33  # set dimensions  
34  L0=1.;L1=1.  if not HAVE_DUDLEY:
35  # bottom temperature:      print("Dudley module not available")
36  T_bot=100  else:
37  # location, size and value of heat source      # set dimensions
38  xc=[0.3,0.4]; r=0.1; Qc=3000      L0=1.;L1=1.
39  # create domain      # bottom temperature:
40  mydomain=Rectangle(l0=L0,l1=L1,n0=20,n1=20)      T_bot=100
41  x=mydomain.getX()      # location, size and value of heat source
42  k=1      xc=[0.3,0.4]; r=0.1; Qc=3000
43  # temperature for boundary condition      # create domain
44  T_D=T_bot/L1*(L1-x[1])      mydomain=Rectangle(l0=L0,l1=L1,n0=20,n1=20)
45  # heat source      x=mydomain.getX()
46  Q=Qc*whereNegative(length(x-xc)-r)      k=1
47  # create PDE:      # temperature for boundary condition
48  mypde=LinearPDE(mydomain)      T_D=T_bot/L1*(L1-x[1])
49  mypde.setSymmetryOn()      # heat source
50  # set coefficients:      Q=Qc*whereNegative(length(x-xc)-r)
51  mypde.setValue(A=k*kronecker(mydomain),Y=Q, r=T_D, \      # create PDE:
52                  q=whereZero(x[1])+whereZero(x[1]-L1))      mypde=LinearPDE(mydomain)
53  # get temperature:      mypde.setSymmetryOn()
54  T=mypde.getSolution()      # set coefficients:
55  # write to file:      mypde.setValue(A=k*kronecker(mydomain),Y=Q, r=T_D, \
56  saveVTK("u.vtu",T=T)                      q=whereZero(x[1])+whereZero(x[1]-L1))
57        # get temperature:
58        T=mypde.getSolution()
59        # write to file:
60        saveVTK("u.vtu",T=T)
61    

Legend:
Removed from v.5287  
changed lines
  Added in v.5288

  ViewVC Help
Powered by ViewVC 1.1.26