/[escript]/trunk/doc/examples/cookbook/example10b.py
ViewVC logotype

Diff of /trunk/doc/examples/cookbook/example10b.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 35  Author: Antony Hallam antony.hallam@uqco Line 35  Author: Antony Hallam antony.hallam@uqco
35  from esys.escript import * # This imports everything from the escript library  from esys.escript import * # This imports everything from the escript library
36  from esys.escript.unitsSI import *  from esys.escript.unitsSI import *
37  from esys.escript.linearPDEs import LinearPDE # This defines LinearPDE as LinearPDE  from esys.escript.linearPDEs import LinearPDE # This defines LinearPDE as LinearPDE
 from esys.finley import Brick # This imports the rectangle domain function from finley  
38  from esys.weipa import saveVTK # This imports the VTK file saver from weipa  from esys.weipa import saveVTK # This imports the VTK file saver from weipa
39  import os, sys #This package is necessary to handle saving our data.  import os, sys #This package is necessary to handle saving our data.
40    
41  from esys.escript.pdetools import Projector, Locator  from esys.escript.pdetools import Projector, Locator
42  import numpy as np  import numpy as np
43    
44    try:
45        # This imports the rectangle domain function
46        from esys.finley import Brick
47        HAVE_FINLEY = True
48    except ImportError:
49        print("Finley module not available")
50        HAVE_FINLEY = False
51  ########################################################MPI WORLD CHECK  ########################################################MPI WORLD CHECK
52  if getMPISizeWorld() > 1:  if getMPISizeWorld() > 1:
53      import sys      import sys
54      print("This example will not run in an MPI world.")      print("This example will not run in an MPI world.")
55      sys.exit(0)      sys.exit(0)
56    
57  #################################################ESTABLISHING VARIABLES  if HAVE_FINLEY:
58  #Domain related.      #################################################ESTABLISHING VARIABLES
59  mx = 500*m #meters - model length      #Domain related.
60  my = 500*m #meters - model width      mx = 500*m #meters - model length
61  mz = -4000*m      my = 500*m #meters - model width
62  ndx = 15 # mesh steps in x direction      mz = -4000*m
63  ndy = 15 # mesh steps in y direction - one dimension means one element      ndx = 15 # mesh steps in x direction
64  ndz = 10      ndy = 15 # mesh steps in y direction - one dimension means one element
65  #PDE related      ndz = 10
66  rho=100.0      #PDE related
67  rholoc=[0,0,-0]      rho=100.0
68  G=6.67300*10E-11      rholoc=[0,0,-0]
69        G=6.67300*10E-11
70  ################################################ESTABLISHING PARAMETERS  
71  #the folder to put our outputs in, leave blank "" for script path      ################################################ESTABLISHING PARAMETERS
72  save_path= os.path.join("data","example10")      #the folder to put our outputs in, leave blank "" for script path
73  #ensure the dir exists      save_path= os.path.join("data","example10")
74  mkDir(save_path)      #ensure the dir exists
75        mkDir(save_path)
76  ####################################################DOMAIN CONSTRUCTION  
77  domain = Brick(l0=mx,l1=my,n0=ndx, n1=ndy,l2=mz,n2=ndz)      ####################################################DOMAIN CONSTRUCTION
78  x=Solution(domain).getX()      domain = Brick(l0=mx,l1=my,n0=ndx, n1=ndy,l2=mz,n2=ndz)
79  x=x-[mx/2,my/2,mz/2]      x=Solution(domain).getX()
80  domain.setX(interpolate(x, ContinuousFunction(domain)))      x=x-[mx/2,my/2,mz/2]
81  mask=wherePositive(100-length(x-rholoc))      domain.setX(interpolate(x, ContinuousFunction(domain)))
82  rho=rho*mask      mask=wherePositive(100-length(x-rholoc))
83  kro=kronecker(domain)      rho=rho*mask
84        kro=kronecker(domain)
85  mass=rho*vol(domain)  
86  ipot=FunctionOnBoundary(domain)      mass=rho*vol(domain)
87  xb=ipot.getX()      ipot=FunctionOnBoundary(domain)
88        xb=ipot.getX()
89  q=whereZero(x[2]-inf(x[2]))  
90  ###############################################ESCRIPT PDE CONSTRUCTION      q=whereZero(x[2]-inf(x[2]))
91        ###############################################ESCRIPT PDE CONSTRUCTION
92  mypde=LinearPDE(domain)  
93  mypde.setValue(A=kro,Y=4.*3.1415*G*rho,q=q,r=0)      mypde=LinearPDE(domain)
94  mypde.setSymmetryOn()      mypde.setValue(A=kro,Y=4.*3.1415*G*rho,q=q,r=0)
95  sol=mypde.getSolution()      mypde.setSymmetryOn()
96  saveVTK(os.path.join(save_path,"ex10b.vtu"),\      sol=mypde.getSolution()
97          grav_pot=sol,\      saveVTK(os.path.join(save_path,"ex10b.vtu"),\
98          g_field=-grad(sol),\              grav_pot=sol,\
99          g_fieldz=-grad(sol)*[0,0,1],\              g_field=-grad(sol),\
100          gz=length(-grad(sol)*[0,0,1]))              g_fieldz=-grad(sol)*[0,0,1],\
101                gz=length(-grad(sol)*[0,0,1]))
102  ################################################MODEL SIZE SAMPLING  
103  sampler=[]      ################################################MODEL SIZE SAMPLING
104  for i in range(-250,250,1):      sampler=[]
105      sampler.append([i,0,250])      for i in range(-250,250,1):
106            sampler.append([i,0,250])
107  sample=[] # array to hold values  
108  rec=Locator(domain,sampler) #location to record      sample=[] # array to hold values
109  psol=rec.getValue(sol)      rec=Locator(domain,sampler) #location to record
110  np.savetxt(os.path.join(save_path,"example10b_%04d.asc"%mx),psol)      psol=rec.getValue(sol)
111        np.savetxt(os.path.join(save_path,"example10b_%04d.asc"%mx),psol)
112    

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

  ViewVC Help
Powered by ViewVC 1.1.26