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

Contents of /trunk/doc/examples/cookbook/heatrefraction_mesher001.py

Parent Directory Parent Directory | Revision Log Revision Log


Revision 2667 - (show annotations)
Thu Sep 17 01:49:11 2009 UTC (9 years, 10 months ago) by jfenwick
File MIME type: text/x-python
File size: 4370 byte(s)
Renamed the main cookbook tex file to match our convention.
Replaced doc/cookbook/figures/heatrefraction002contqu.pdf with
a version which is actually pdf. However it needs to be regnerated since
it it sideways.

The examples have had their copyright notices fixed (dates were too early).
sb2.py has been removed since it uses pyvisi.

scons will now build the cookbook as parts of a docs build.
Also in reposnse to :
scons cookbook_pdf


1
2 ########################################################
3 #
4 # Copyright (c) 2009 by University of Queensland
5 # Earth Systems Science Computational Center (ESSCC)
6 # http://www.uq.edu.au/esscc
7 #
8 # Primary Business: Queensland, Australia
9 # Licensed under the Open Software License version 3.0
10 # http://www.opensource.org/licenses/osl-3.0.php
11 #
12 ########################################################
13
14 __copyright__="""Copyright (c) 2009 by University of Queensland
15 Earth Systems Science Computational Center (ESSCC)
16 http://www.uq.edu.au/esscc
17 Primary Business: Queensland, Australia"""
18 __license__="""Licensed under the Open Software License version 3.0
19 http://www.opensource.org/licenses/osl-3.0.php"""
20 __url__="https://launchpad.net/escript-finley"
21
22 """
23 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
33 from esys.pycad.gmsh import Design #Finite Element meshing package
34 from esys.finley import MakeDomain #Converter for escript
35 import os #file path tool
36 import numpy as np #numerial python for arrays
37 from math import * # math package
38
39 #used to construct polygons for plotting from pycad
40 from cblib import getLoopCoords, needdirs
41
42
43
44 #################################################ESTABLISHING VARIABLES
45 #set modal to 1 for a syncline or -1 for an anticline structural
46 #configuration
47 modal=1
48
49 # the folder to put our outputs in, leave blank "" for script path -
50 # note this folder path must exist to work
51 save_path= os.path.join("data","heatrefrac001")
52 needdirs([save_path])
53
54 ################################################ESTABLISHING PARAMETERS
55 #Model Parameters
56 width=5000.0 #width of model
57 depth=-6000.0 #depth of model
58 sspl=51 #number of discrete points in spline
59 dsp=width/(sspl-1) #dx of spline steps for width
60 dep_sp=2500.0 #avg depth of spline
61 h_sp=1500.0 #heigh of spline
62 orit=-1.0 #orientation of spline 1.0=>up -1.0=>down
63
64 ####################################################DOMAIN CONSTRUCTION
65 # Domain Corners
66 p0=Point(0.0, 0.0, 0.0)
67 p1=Point(0.0, depth, 0.0)
68 p2=Point(width, depth, 0.0)
69 p3=Point(width, 0.0, 0.0)
70 # Join corners in anti-clockwise manner.
71 l01=Line(p0, p1)
72 l12=Line(p1, p2)
73 l23=Line(p2, p3)
74 l30=Line(p3, p0)
75 # Join line segments to create domain boundary.
76 c=CurveLoop(l01, l12, l23, l30)
77 # Generate Material Boundary
78 x=[ Point(i*dsp\
79 ,-dep_sp+modal*orit*h_sp*cos(pi*i*dsp/dep_sp+pi))\
80 for i in range(0,sspl)\
81 ]
82 mysp = Spline(*tuple(x))
83 # Start and end of material boundary.
84 x1=Spline.getStartPoint(mysp)
85 x2=Spline.getEndPoint(mysp)
86
87 # Create TOP BLOCK
88 # lines
89 tbl1=Line(p0,x1)
90 tbl2=mysp
91 tbl3=Line(x2,p3)
92 # curve
93 tblockloop = CurveLoop(tbl1,tbl2,tbl3,l30)
94 # surface
95 tblock = PlaneSurface(tblockloop)
96 # polygon for python plotting
97 tpg = getLoopCoords(tblockloop)
98 np.savetxt(os.path.join(save_path,"toppg"),tpg,delimiter=" ")
99
100 # Create BOTTOM BLOCK
101 # lines
102 bbl1=Line(x1,p1)
103 bbl3=Line(p2,x2)
104 bbl4=-mysp
105 # curve
106 bblockloop = CurveLoop(bbl1,l12,bbl3,bbl4)
107 # surface
108 bblock = PlaneSurface(bblockloop)
109
110 #clockwise check as splines must be set as polygons in the point order
111 #they were created. Otherwise get a line across plot.
112 bblockloop2=CurveLoop(mysp,Line(x2,p2),Line(p2,p1),Line(p1,x1))
113 bpg = getLoopCoords(bblockloop2)
114 np.savetxt(os.path.join(save_path,"botpg"),bpg,delimiter=" ")
115
116 #############################################EXPORTING MESH FOR ESCRIPT
117 # Create a Design which can make the mesh
118 d=Design(dim=2, element_size=200)
119 # Add the subdomains and flux boundaries.
120 d.addItems(PropertySet("top",tblock),PropertySet("bottom",bblock),\
121 PropertySet("linebottom",l12))
122 # Create the geometry, mesh and Escript domain
123 d.setScriptFileName(os.path.join(save_path,\
124 "heatrefraction_mesh001.geo"))
125
126 d.setMeshFileName(os.path.join(save_path,"heatrefraction_mesh001.msh"))
127 domain=MakeDomain(d, integrationOrder=-1, reducedIntegrationOrder=-1,\
128 optimizeLabeling=True)
129 # Create a file that can be read back in to python with
130 # mesh=ReadMesh(fileName)
131 domain.write(os.path.join(save_path,"heatrefraction_mesh001.fly"))
132
133

  ViewVC Help
Powered by ViewVC 1.1.26