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

Diff of /trunk/doc/examples/cookbook/heatrefraction_mesher002.py

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

revision 2657 by jfenwick, Mon Sep 7 00:06:15 2009 UTC revision 2658 by ahallam, Thu Sep 10 02:58:44 2009 UTC
# Line 22  __url__="https://launchpad.net/escript-f Line 22  __url__="https://launchpad.net/escript-f
22  """  """
23  Author: Antony Hallam antony.hallam@uqconnect.edu.au  Author: Antony Hallam antony.hallam@uqconnect.edu.au
24  """  """
25    ############################################################FILE HEADER
26    # heatrefraction_mesher002.py
27    # Create a 3 block fault and overburden style model in 2d using pycad
28    # meshing tools.
29    
30    #######################################################EXTERNAL MODULES
31  from esys.pycad import *  from esys.pycad import *
32  from esys.pycad.gmsh import Design  from esys.pycad.gmsh import Design
33  from esys.finley import MakeDomain  from esys.finley import MakeDomain
34  import os  import os
35  import numpy as np  import numpy as np
36  from math import *  from math import *
37  from cblib import needdirs  from cblib import needdirs, getLoopCoords
38    
39    #################################################ESTABLISHING VARIABLES
40  # routine to find consecutive coordinates of a loop in pycad  # where to put output files
41  def loopcoords(loop):  save_path= os.path.join("data","heatrefrac002")
     # return all construction points of input  
     temp = loop.getConstructionPoints()  
     #create a numpy array for xyz components or construction points  
     coords = np.zeros([len(temp),3],float)  
     #place construction points in array  
     for i in range(0,len(temp)):  
         coords[i,:]=temp[i].getCoordinates()  
     #return a numpy array  
     return coords  
   
 save_path = "data/heatrefrac002"  
42  needdirs([save_path])  needdirs([save_path])
43    
44    ################################################ESTABLISHING PARAMETERS
45  # Overall Domain  # Overall Domain
46  p0=Point(0.0,        0.0, 0.0)  p0=Point(0.0,        0.0, 0.0)
47  p1=Point(0.0,    -6000.0, 0.0)  p1=Point(0.0,    -6000.0, 0.0)
48  p2=Point(5000.0, -6000.0, 0.0)  p2=Point(5000.0, -6000.0, 0.0)
49  p3=Point(5000.0,     0.0, 0.0)  p3=Point(5000.0,     0.0, 0.0)
50    
51    ####################################################DOMAIN CONSTRUCTION
52  l01=Line(p0, p1)  l01=Line(p0, p1)
53  l12=Line(p1, p2)  l12=Line(p1, p2)
54  l23=Line(p2, p3)  l23=Line(p2, p3)
# Line 60  l30=Line(p3, p0) Line 56  l30=Line(p3, p0)
56    
57  c=CurveLoop(l01, l12, l23, l30)  c=CurveLoop(l01, l12, l23, l30)
58    
59  # Material Boundary  # Generate Material Boundary
   
60  p4=Point(0.0,    -2400.0, 0.0)  p4=Point(0.0,    -2400.0, 0.0)
61  p5=Point(2000.0, -2400.0, 0.0)  p5=Point(2000.0, -2400.0, 0.0)
62  p6=Point(3000.0, -6000.0, 0.0)  p6=Point(3000.0, -6000.0, 0.0)
63  p7=Point(5000.0, -2400.0, 0.0)  p7=Point(5000.0, -2400.0, 0.0)
64            
65  # TOP BLOCK  # Create TOP BLOCK
66  tbl1=Line(p0,p4)  tbl1=Line(p0,p4)
67  tbl2=Line(p4,p5)  tbl2=Line(p4,p5)
68  tbl3=Line(p5,p7)  tbl3=Line(p5,p7)
69  tbl4=Line(p7,p3)  tbl4=Line(p7,p3)
70  tblockloop = CurveLoop(tbl1,tbl2,tbl3,tbl4,l30)  tblockloop = CurveLoop(tbl1,tbl2,tbl3,tbl4,l30)
71  tblock = PlaneSurface(tblockloop)  tblock = PlaneSurface(tblockloop)
72    # python plotting polygon
73  tpg = loopcoords(tblockloop)  tpg = loopcoords(tblockloop)
74  np.savetxt(os.path.join(save_path,"toppg"),tpg,delimiter=" ")  np.savetxt(os.path.join(save_path,"toppg"),tpg,delimiter=" ")
75    
76  # BOTTOM BLOCK LEFT  # Create BOTTOM BLOCK LEFT
77  bbll1=Line(p4,p1)  bbll1=Line(p4,p1)
78  bbll2=Line(p1,p6)  bbll2=Line(p1,p6)
79  bbll3=Line(p6,p5)  bbll3=Line(p6,p5)
# Line 91  bblockl = PlaneSurface(bblockloopl) Line 86  bblockl = PlaneSurface(bblockloopl)
86  bpg = loopcoords(bblockloopl)  bpg = loopcoords(bblockloopl)
87  np.savetxt(os.path.join(save_path,"botpgl"),bpg,delimiter=" ")  np.savetxt(os.path.join(save_path,"botpgl"),bpg,delimiter=" ")
88    
89  # BOTTOM BLOCK RIGHT  # Create BOTTOM BLOCK RIGHT
90  bbrl1=Line(p6,p2)  bbrl1=Line(p6,p2)
91  bbrl2=Line(p2,p7)  bbrl2=Line(p2,p7)
92  bbrl3=-tbl3  bbrl3=-tbl3
# Line 104  bblockr = PlaneSurface(bblockloopr) Line 99  bblockr = PlaneSurface(bblockloopr)
99  bpg = loopcoords(bblockloopr)  bpg = loopcoords(bblockloopr)
100  np.savetxt(os.path.join(save_path,"botpgr"),bpg,delimiter=" ")  np.savetxt(os.path.join(save_path,"botpgr"),bpg,delimiter=" ")
101    
102    #############################################EXPORTING MESH FOR ESCRIPT
103  # Create a Design which can make the mesh  # Create a Design which can make the mesh
104  d=Design(dim=2, element_size=200)  d=Design(dim=2, element_size=200)
105  # Add the subdomains and flux boundaries.  # Add the subdomains and flux boundaries.
106  d.addItems(PropertySet("top",tblock),PropertySet("bottomleft",bblockl),PropertySet("bottomright",bblockr),PropertySet("linebottom",bbll21, bbrl1))  d.addItems(PropertySet("top",tblock),\
107                           PropertySet("bottomleft",bblockl),\
108                           PropertySet("bottomright",bblockr),\
109                           PropertySet("linebottom",bbll21, bbrl1))
110    
111  # Create the geometry, mesh and Escript domain  # Create the geometry, mesh and Escript domain
112  d.setScriptFileName(os.path.join(save_path,"heatrefraction_mesh003.geo"))  d.setScriptFileName(os.path.join(save_path,\
113                          "heatrefraction_mesh003.geo"))
114    
115  d.setMeshFileName(os.path.join(save_path,"heatrefraction_mesh003.msh"))  d.setMeshFileName(os.path.join(save_path,"heatrefraction_mesh003.msh"))
116  domain=MakeDomain(d, integrationOrder=-1, reducedIntegrationOrder=-1, optimizeLabeling=True)  domain=MakeDomain(d, integrationOrder=-1,\
117  # Create a file that can be read back in to python with mesh=ReadMesh(fileName)                       reducedIntegrationOrder=-1,\
118                         optimizeLabeling=True)
119    # Create a file that can be read back in to python with
120    # mesh=ReadMesh(fileName)
121  domain.write(os.path.join(save_path,"heatrefraction_mesh003.fly"))  domain.write(os.path.join(save_path,"heatrefraction_mesh003.fly"))
122    

Legend:
Removed from v.2657  
changed lines
  Added in v.2658

  ViewVC Help
Powered by ViewVC 1.1.26