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

Diff of /trunk/doc/examples/cookbook/heatrefraction_mesher001.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 23  __url__="https://launchpad.net/escript-f Line 23  __url__="https://launchpad.net/escript-f
23  Author: Antony Hallam antony.hallam@uqconnect.edu.au  Author: Antony Hallam antony.hallam@uqconnect.edu.au
24  """  """
25    
26    ############################################################FILE HEADER
27    # heatrefraction_mesher001.py
28    # Create either a 2D syncline or anticline model using pycad meshing
29    # tools.
30    
31    #######################################################EXTERNAL MODULES
32  from esys.pycad import * #domain constructor  from esys.pycad import * #domain constructor
33  from esys.pycad.gmsh import Design #Finite Element meshing package  from esys.pycad.gmsh import Design #Finite Element meshing package
34  from esys.finley import MakeDomain #Converter for escript  from esys.finley import MakeDomain #Converter for escript
# Line 32  from math import * # math package Line 38  from math import * # math package
38  #used to construct polygons for plotting from pycad  #used to construct polygons for plotting from pycad
39  from cblib import getLoopCoords, needdirs  from cblib import getLoopCoords, needdirs
40    
41    #################################################ESTABLISHING VARIABLES
42  #set modal to 1 for a syncline or -1 for an anticline structural  #set modal to 1 for a syncline or -1 for an anticline structural
43  #configuration    #configuration  
44  modal=1  modal=1
45    
46  # the folder to put our outputs in, leave blank "" for script path -  # the folder to put our outputs in, leave blank "" for script path -
47  # note this folder path must exist to work  # note this folder path must exist to work
48  save_path = "data/heatrefrac001"  save_path= os.path.join("data","heatrefrac001")
49  needdirs([save_path])  needdirs([save_path])
50    
51    ################################################ESTABLISHING PARAMETERS
52  #Model Parameters  #Model Parameters
53  width=5000.0   #width of model  width=5000.0   #width of model
54  depth=-6000.0  #depth of model  depth=-6000.0  #depth of model
# Line 50  dep_sp=2500.0 #avg depth of spline Line 58  dep_sp=2500.0 #avg depth of spline
58  h_sp=1500.0 #heigh of spline  h_sp=1500.0 #heigh of spline
59  orit=-1.0 #orientation of spline 1.0=>up -1.0=>down  orit=-1.0 #orientation of spline 1.0=>up -1.0=>down
60    
61  # Overall Domain  ####################################################DOMAIN CONSTRUCTION
62    # Domain Corners
63  p0=Point(0.0,      0.0, 0.0)  p0=Point(0.0,      0.0, 0.0)
64  p1=Point(0.0,    depth, 0.0)  p1=Point(0.0,    depth, 0.0)
65  p2=Point(width, depth, 0.0)  p2=Point(width, depth, 0.0)
66  p3=Point(width,   0.0, 0.0)  p3=Point(width,   0.0, 0.0)
67    # Join corners in anti-clockwise manner.
68  l01=Line(p0, p1)  l01=Line(p0, p1)
69  l12=Line(p1, p2)  l12=Line(p1, p2)
70  l23=Line(p2, p3)  l23=Line(p2, p3)
71  l30=Line(p3, p0)  l30=Line(p3, p0)
72    # Join line segments to create domain boundary.
73  c=CurveLoop(l01, l12, l23, l30)  c=CurveLoop(l01, l12, l23, l30)
74    # Generate Material Boundary
 # Material Boundary  
75  x=[ Point(i*dsp\  x=[ Point(i*dsp\
76      ,-dep_sp+modal*orit*h_sp*cos(pi*i*dsp/dep_sp+pi))\      ,-dep_sp+modal*orit*h_sp*cos(pi*i*dsp/dep_sp+pi))\
77       for i in range(0,sspl)\       for i in range(0,sspl)\
78    ]    ]
       
79  mysp = Spline(*tuple(x))  mysp = Spline(*tuple(x))
80    # Start and end of material boundary.
81  x1=Spline.getStartPoint(mysp)  x1=Spline.getStartPoint(mysp)
82  x2=Spline.getEndPoint(mysp)  x2=Spline.getEndPoint(mysp)
83            
84  # TOP BLOCK  #  Create TOP BLOCK
85    # lines
86  tbl1=Line(p0,x1)  tbl1=Line(p0,x1)
87  tbl2=mysp  tbl2=mysp
88  tbl3=Line(x2,p3)  tbl3=Line(x2,p3)
89    # curve
90  tblockloop = CurveLoop(tbl1,tbl2,tbl3,l30)  tblockloop = CurveLoop(tbl1,tbl2,tbl3,l30)
91    # surface
92  tblock = PlaneSurface(tblockloop)  tblock = PlaneSurface(tblockloop)
93    # polygon for python plotting
94  tpg = getLoopCoords(tblockloop)  tpg = getLoopCoords(tblockloop)
95  np.savetxt(os.path.join(save_path,"toppg"),tpg,delimiter=" ")  np.savetxt(os.path.join(save_path,"toppg"),tpg,delimiter=" ")
96    
97  # BOTTOM BLOCK  # Create BOTTOM BLOCK
98    # lines
99  bbl1=Line(x1,p1)  bbl1=Line(x1,p1)
100  bbl3=Line(p2,x2)  bbl3=Line(p2,x2)
101  bbl4=-mysp  bbl4=-mysp
102    # curve
103  bblockloop = CurveLoop(bbl1,l12,bbl3,bbl4)  bblockloop = CurveLoop(bbl1,l12,bbl3,bbl4)
104    # surface
105  bblock = PlaneSurface(bblockloop)  bblock = PlaneSurface(bblockloop)
106    
107  #clockwise check as splines must be set as polygons in the point order  #clockwise check as splines must be set as polygons in the point order
# Line 96  bblockloop2=CurveLoop(mysp,Line(x2,p2),L Line 110  bblockloop2=CurveLoop(mysp,Line(x2,p2),L
110  bpg = getLoopCoords(bblockloop2)  bpg = getLoopCoords(bblockloop2)
111  np.savetxt(os.path.join(save_path,"botpg"),bpg,delimiter=" ")  np.savetxt(os.path.join(save_path,"botpg"),bpg,delimiter=" ")
112    
113    #############################################EXPORTING MESH FOR ESCRIPT
114  # Create a Design which can make the mesh  # Create a Design which can make the mesh
115  d=Design(dim=2, element_size=200)  d=Design(dim=2, element_size=200)
116  # Add the subdomains and flux boundaries.  # Add the subdomains and flux boundaries.
117  d.addItems(PropertySet("top",tblock),PropertySet("bottom",bblock),PropertySet("linebottom",l12))  d.addItems(PropertySet("top",tblock),PropertySet("bottom",bblock),\
118                                         PropertySet("linebottom",l12))
119  # Create the geometry, mesh and Escript domain  # Create the geometry, mesh and Escript domain
120  d.setScriptFileName(os.path.join(save_path,"heatrefraction_mesh001.geo"))  d.setScriptFileName(os.path.join(save_path,\
121                                       "heatrefraction_mesh001.geo"))
122    
123  d.setMeshFileName(os.path.join(save_path,"heatrefraction_mesh001.msh"))  d.setMeshFileName(os.path.join(save_path,"heatrefraction_mesh001.msh"))
124  domain=MakeDomain(d, integrationOrder=-1, reducedIntegrationOrder=-1, optimizeLabeling=True)  domain=MakeDomain(d, integrationOrder=-1, reducedIntegrationOrder=-1,\
125  # Create a file that can be read back in to python with mesh=ReadMesh(fileName)                     optimizeLabeling=True)
126    # Create a file that can be read back in to python with
127    # mesh=ReadMesh(fileName)
128  domain.write(os.path.join(save_path,"heatrefraction_mesh001.fly"))  domain.write(os.path.join(save_path,"heatrefraction_mesh001.fly"))
129    
130    

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

  ViewVC Help
Powered by ViewVC 1.1.26