# Diff of /trunk/pyvisi/test/python/gen_data.py

revision 850 by gross, Mon Aug 14 05:36:17 2006 UTC revision 851 by gross, Mon Sep 18 07:28:18 2006 UTC
11  from esys.escript import *  from esys.escript import *
12    from esys.escript.pdetools import Projector
13  from esys import finley  from esys import finley
14  #  #
15  #  #
# Line 23  NE2=int(NE0/L0*L2) Line 24  NE2=int(NE0/L0*L2)
24  #  #
25  #  generate domains:  #  generate domains:
26  #  #
27  if DIM == 2:  for DIM in [2,3]:
28   domain=finley.Rectangle(NE0,NE1,1,l0=L0,l1=L1)    if DIM == 2:
29   trafo=numarray.array([[0.,1.],[-1.,0.]])     domain=finley.Rectangle(NE0,NE1,1,l0=L0,l1=L1)
30   cen=numarray.array([L0/2.,L1/2.])     trafo=numarray.array([[0.,1.],[-1.,0.]])
31  else:     cen=numarray.array([L0/2.,L1/2.])
32   domain=finley.Brick(NE0,NE1,NE2,1,l0=L0,l1=L1,l2=L2)    else:
33   trafo=numarray.array([[0.,1.,0.],[-1.,0.,0.],[0.,0.,1.]])     domain=finley.Brick(NE0,NE1,NE2,1,l0=L0,l1=L1,l2=L2)
34   cen=numarray.array([L0/2.,L1/2.,L2/2.])     trafo=numarray.array([[0.,1.,0.],[-1.,0.,0.],[0.,0.,1.]])
35  #     cen=numarray.array([L0/2.,L1/2.,L2/2.])
36  # get function spaces:    pp=Projector(domain)
37  #    #
38  c=ContinuousFunction(domain)    # get function spaces:
39  f=Function(domain)    #
40  b=FunctionOnBoundary(domain)    c=ContinuousFunction(domain)
41  #    f=Function(domain)
42  # get coordinates    b=FunctionOnBoundary(domain)
43  #    #
44  c_x=c.getX()-cen    # get coordinates
45  c_r=length(c_x)    #
46      c_x=c.getX()-cen
47  f_x=f.getX()-cen    c_r=length(c_x)
48  f_r=length(f_x)    c_t=matrix_mult(trafo,c_x/(c_r+1.e-15))
49
50  b_x=b.getX()-cen    f_x=f.getX()-cen
51  b_r=length(b_x)    f_r=length(f_x)
52  #
53  #      b_x=b.getX()-cen
54  #    b_r=length(b_x)
55  saveVTK("interior_%dD.xml"%DIM,temperature=sin(2*c_r),temperature_cell=sin(2*f_r), \    #
56                                 velocity=matrix_mult(trafo,c_x/(c_r+1.e-15)), velocity_cell=matrix_mult(trafo,f_x/f_r))    #
57  saveVTK("boundary_%dD.xml"%DIM,temperature=sin(b_r),velocity=matrix_mult(trafo,b_x/(b_r+1.e-15)))    #
58      saveVTK("interior_%dD.xml"%DIM,temperature=sin(2*c_r),temperature_cell=sin(2*f_r), \
59                                     velocity=matrix_mult(trafo,c_x/(c_r+1.e-15)), velocity_cell=matrix_mult(trafo,f_x/f_r), \