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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 2914 - (show annotations)
Wed Feb 3 05:38:12 2010 UTC (10 years ago) by gross
File MIME type: text/x-python
File size: 3922 byte(s)
missing import added.
1
2 ########################################################
3 #
4 # Copyright (c) 2009-2010 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-2010 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 ############################################################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 *
32 from esys.pycad.gmsh import Design
33 from esys.finley import MakeDomain
34 import os
35 import numpy as np
36 from math import *
37 from cblib1 import getLoopCoords
38 from esys.escript import mkDir
39
40 #################################################ESTABLISHING VARIABLES
41 # where to put output files
42 save_path= os.path.join("data","heatrefrac002")
43 mkDir(save_path)
44
45 ################################################ESTABLISHING PARAMETERS
46 # Overall Domain
47 p0=Point(0.0, 0.0, 0.0)
48 p1=Point(0.0, -6000.0, 0.0)
49 p2=Point(5000.0, -6000.0, 0.0)
50 p3=Point(5000.0, 0.0, 0.0)
51
52 ####################################################DOMAIN CONSTRUCTION
53 l01=Line(p0, p1)
54 l12=Line(p1, p2)
55 l23=Line(p2, p3)
56 l30=Line(p3, p0)
57
58 c=CurveLoop(l01, l12, l23, l30)
59
60 # Generate Material Boundary
61 p4=Point(0.0, -2400.0, 0.0)
62 p5=Point(2000.0, -2400.0, 0.0)
63 p6=Point(3000.0, -6000.0, 0.0)
64 p7=Point(5000.0, -2400.0, 0.0)
65
66 # Create TOP BLOCK
67 tbl1=Line(p0,p4)
68 tbl2=Line(p4,p5)
69 tbl3=Line(p5,p7)
70 tbl4=Line(p7,p3)
71 tblockloop = CurveLoop(tbl1,tbl2,tbl3,tbl4,l30)
72 tblock = PlaneSurface(tblockloop)
73 # python plotting polygon
74 tpg = getLoopCoords(tblockloop)
75 np.savetxt(os.path.join(save_path,"toppg"),tpg,delimiter=" ")
76
77 # Create BOTTOM BLOCK LEFT
78 bbll1=Line(p4,p1)
79 bbll2=Line(p1,p6)
80 bbll3=Line(p6,p5)
81 bbll4=-tbl2
82 bblockloopl = CurveLoop(bbll1,bbll2,bbll3,bbll4)
83 bblockl = PlaneSurface(bblockloopl)
84
85 #clockwise check
86 #bblockloopl2=CurveLoop(mysp,Line(x2,p2),Line(p2,p1),Line(p1,x1))
87 bpg = getLoopCoords(bblockloopl)
88 np.savetxt(os.path.join(save_path,"botpgl"),bpg,delimiter=" ")
89
90 # Create BOTTOM BLOCK RIGHT
91 bbrl1=Line(p6,p2)
92 bbrl2=Line(p2,p7)
93 bbrl3=-tbl3
94 bbrl4=-bbll3
95 bblockloopr = CurveLoop(bbrl1,bbrl2,bbrl3,bbrl4)
96 bblockr = PlaneSurface(bblockloopr)
97
98 #clockwise check
99 #bblockloopr2=CurveLoop(mysp,Line(x2,p2),Line(p2,p1),Line(p1,x1))
100 bpg = getLoopCoords(bblockloopr)
101 np.savetxt(os.path.join(save_path,"botpgr"),bpg,delimiter=" ")
102
103 #############################################EXPORTING MESH FOR ESCRIPT
104 # Create a Design which can make the mesh
105 d=Design(dim=2, element_size=200)
106 # Add the subdomains and flux boundaries.
107 d.addItems(PropertySet("top",tblock),\
108 PropertySet("bottomleft",bblockl),\
109 PropertySet("bottomright",bblockr),\
110 PropertySet("linebottom",bbll2, bbrl1))
111
112 # Create the geometry, mesh and Escript domain
113 d.setScriptFileName(os.path.join(save_path,\
114 "heatrefraction_mesh003.geo"))
115
116 d.setMeshFileName(os.path.join(save_path,"heatrefraction_mesh003.msh"))
117 domain=MakeDomain(d, integrationOrder=-1,\
118 reducedIntegrationOrder=-1,\
119 optimizeLabeling=True)
120 # Create a file that can be read back in to python with
121 # mesh=ReadMesh(fileName)
122 domain.write(os.path.join(save_path,"heatrefraction_mesh003.fly"))
123

  ViewVC Help
Powered by ViewVC 1.1.26