/[escript]/trunk/pyvisi/test/python/run_escript_with_lazy_evaluation.py
ViewVC logotype

Diff of /trunk/pyvisi/test/python/run_escript_with_lazy_evaluation.py

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

revision 1165 by jongui, Wed May 23 04:17:52 2007 UTC revision 1166 by jongui, Thu May 24 00:56:00 2007 UTC
# Line 75  class TestEscriptMap(unittest.TestCase, Line 75  class TestEscriptMap(unittest.TestCase,
75          dc = DataCollector(source = Source.ESCRIPT)          dc = DataCollector(source = Source.ESCRIPT)
76    
77          # Create a Map.          # Create a Map.
78          m = Map(scene = s, data_collector = dc,          m = Map(scene = s, data_collector = dc, \
79                  viewport = Viewport.SOUTH_WEST, lut = Lut.COLOR,                  viewport = Viewport.SOUTH_WEST, lut = Lut.COLOR, \
80                  cell_to_point = False, outline = True)                  cell_to_point = False, outline = True)
81    
82          # Create a Camera.          # Create a Camera.
83          c = Camera(scene = s, viewport = Viewport.SOUTH_WEST)          c = Camera(scene = s, viewport = Viewport.SOUTH_WEST)
84    
85          # ... start iteration:          # ... start iteration:
86          while t<0.2:          while t<0.4:
87              i+=1              i+=1
88              t+=h              t+=h
89              mypde.setValue(Y=qH+rhocp/h*T)              mypde.setValue(Y=qH+rhocp/h*T)
# Line 107  class TestEscriptVelocity(unittest.TestC Line 107  class TestEscriptVelocity(unittest.TestC
107          beta=8.          beta=8.
108          T_ref=0.          T_ref=0.
109          T_0=1.          T_0=1.
110    
111          #... generate domain ...          #... generate domain ...
112          mydomain = Brick(l0=1.,l1=1., l2=1.,n0=10, n1=10, n2=10)          mydomain = Brick(l0=1.,l1=1., l2=1.,n0=10, n1=10, n2=10)
113          x=mydomain.getX()          x=mydomain.getX()
# Line 115  class TestEscriptVelocity(unittest.TestC Line 116  class TestEscriptVelocity(unittest.TestC
116          #... open symmetric PDE ...          #... open symmetric PDE ...
117          mypde=LinearPDE(mydomain)          mypde=LinearPDE(mydomain)
118          mypde.setSymmetryOn()          mypde.setSymmetryOn()
119    
120          #... set coefficients ...          #... set coefficients ...
121          C=Tensor4(0.,Function(mydomain))          C=Tensor4(0.,Function(mydomain))
122          for i in range(mydomain.getDim()):          for i in range(mydomain.getDim()):
123            for j in range(mydomain.getDim()):              for j in range(mydomain.getDim()):
124               C[i,i,j,j]+=lam                  C[i,i,j,j]+=lam
125               C[j,i,j,i]+=mu                  C[j,i,j,i]+=mu
126               C[j,i,i,j]+=mu                  C[j,i,i,j]+=mu
127    
128          msk=whereZero(x[0])*[1.,0.,0.] \          msk=whereZero(x[0])*[1.,0.,0.] \
129             +whereZero(x[1])*[0.,1.,0.] \                  +whereZero(x[1])*[0.,1.,0.] \
130             +whereZero(x[2])*[0.,0.,1.]                  +whereZero(x[2])*[0.,0.,1.]
131          sigma0=(lam+2./3.*mu)*alpha*(T-T_ref)*kronecker(mydomain)          sigma0=(lam+2./3.*mu)*alpha*(T-T_ref)*kronecker(mydomain)
132          mypde.setValue(A=C,X=sigma0,q=msk)          mypde.setValue(A=C,X=sigma0,q=msk)
133    
134          #... solve pde ...          #... solve pde ...
135          u=mypde.getSolution()          u=mypde.getSolution()
136          #==============          #... calculate von-Misses
         # ... calculate von-Misses  
         u=mydomain.getX()  
137          g=grad(u)          g=grad(u)
138          sigma=mu*(g+transpose(g))+lam*trace(g)*kronecker(mydomain)-sigma0          sigma=mu*(g+transpose(g))+lam*trace(g)*kronecker(mydomain)-sigma0
139          sigma_mises=\          sigma_mises=sqrt(((sigma[0,0]-sigma[1,1])**2+\
140                  sqrt(((sigma[0,0]-sigma[1,1])**2+(sigma[1,1]-sigma[2,2])**2+ \                  (sigma[1,1]-sigma[2,2])**2+ \
141                  (sigma[2,2]-sigma[0,0])**2)/6. \                  (sigma[2,2]-sigma[0,0])**2)/6. \
142                  +sigma[0,1]**2 + sigma[1,2]**2 + sigma[2,0]**2)                  +sigma[0,1]**2 + sigma[1,2]**2 + sigma[2,0]**2)
143            
144          # Create a Scene.          # Create a Scene.
145          s = Scene(renderer = JPG_RENDERER, x_size = X_SIZE, y_size = Y_SIZE)          s = Scene(renderer = JPG_RENDERER, x_size = X_SIZE, y_size = Y_SIZE)
146          self.scene = s          self.scene = s
147    
148          # Create a DataCollector reading directly from an escript object.          # Create a DataCollector reading directly from an escript object.
149          dc = DataCollector(source = Source.ESCRIPT)          dc = DataCollector(source = Source.ESCRIPT)
         T=length(mydomain.getX())  
150          dc.setData(disp = u, stress = sigma_mises)          dc.setData(disp = u, stress = sigma_mises)
151    
152          # Create a Velocity.          # Create a Velocity.
153          v = Velocity(scene = s, data_collector = dc,          v = Velocity(scene = s, data_collector = dc, \
154                  viewport = Viewport.SOUTH_WEST,                  viewport = Viewport.SOUTH_WEST, \
155                  arrow = Arrow.THREE_D, color_mode = ColorMode.SCALAR,                  arrow = Arrow.THREE_D, color_mode = ColorMode.SCALAR, \
156                  lut = Lut.COLOR, cell_to_point = True, outline = True)                  lut = Lut.COLOR, cell_to_point = True, outline = True)
157          v.setScaleFactor(scale_factor = 0.3)          v.setScaleFactor(scale_factor = 0.3)
158    
# Line 197  class TestEscriptEllipsoid(unittest.Test Line 198  class TestEscriptEllipsoid(unittest.Test
198          mypde.setValue(D=kronecker*rho)          mypde.setValue(D=kronecker*rho)
199          # ... set initial values ....          # ... set initial values ....
200          n=0          n=0
201          # initial value of displacement at point source is constant          # initial value of displacement at point source is constant (U0=0.01)
202          # (U0=0.01) for first two time steps          # for first two time steps
203          u=U0*whereNegative(length(x-xc)-src_radius)*dunit          u=U0*whereNegative(length(x-xc)-src_radius)*dunit
204          u_last=U0*whereNegative(length(x-xc)-src_radius)*dunit          u_last=U0*whereNegative(length(x-xc)-src_radius)*dunit
205          t=0          t=0
# Line 225  class TestEscriptEllipsoid(unittest.Test Line 226  class TestEscriptEllipsoid(unittest.Test
226    
227          # Create an Ellipsoid.          # Create an Ellipsoid.
228          e = Ellipsoid(scene = s, data_collector = dc,          e = Ellipsoid(scene = s, data_collector = dc,
229                 viewport = Viewport.SOUTH_WEST,          viewport = Viewport.SOUTH_WEST,
230                 lut = Lut.COLOR, cell_to_point = True, outline = True)          lut = Lut.COLOR, cell_to_point = True, outline = True)
231          e.setScaleFactor(scale_factor = 0.7)          e.setScaleFactor(scale_factor = 0.7)
232          e.setMaxScaleFactor(max_scale_factor = 100)          e.setMaxScaleFactor(max_scale_factor = 1000)
233          e.setRatio(ratio = 10)          e.setRatio(ratio = 10)
234    
235          # Create a Camera.          # Create a Camera.
236          c = Camera(scene = s, viewport = Viewport.SOUTH_WEST)          c = Camera(scene = s, viewport = Viewport.SOUTH_WEST)
237          c.isometricView()          c.isometricView()
238    
239          while t<0.2:          while t<0.5:
240              # ... get current stress ....              # ... get current stress ....
241              g=grad(u)              g=grad(u)
242              stress=lam*trace(g)*kronecker+mu*(g+transpose(g))              stress=lam*trace(g)*kronecker+mu*(g+transpose(g))
# Line 255  class TestEscriptEllipsoid(unittest.Test Line 256  class TestEscriptEllipsoid(unittest.Test
256              u_pc_y=u_pc[1]              u_pc_y=u_pc[1]
257              u_pc_z=u_pc[2]              u_pc_z=u_pc[2]
258    
259              # ... save current acceleration in units of gravity and displacements              # ... save current acceleration in units of gravity and
260              # === syntetic data              # displacements
             u=domain.getX()  
261              if n==1 or n%10==0:              if n==1 or n%10==0:
262                  dc.setData(acceleration = length(a)/9.81, displacement = u,                  dc.setData(acceleration = length(a)/9.81, displacement = u,
263                       tensor = stress, Ux = u[0])                          tensor = stress, Ux = u[0])
264    
265               # Render the object.                  # Render the object.
266              self.render("TestEscriptEllipsoid%02d.jpg" % (t/10))                  self.render("TestEscriptEllipsoid%02d.jpg" % (n/10))
267    
268    
269  if __name__ == '__main__':  if __name__ == '__main__':

Legend:
Removed from v.1165  
changed lines
  Added in v.1166

  ViewVC Help
Powered by ViewVC 1.1.26