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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 2597 - (show annotations)
Thu Aug 6 03:30:09 2009 UTC (9 years, 8 months ago) by ahallam
File MIME type: text/x-python
File size: 3163 byte(s)
Updates to: cookbook documentation; heat refraction problems; mencoder investigation for *.png animation
1
2 ########################################################
3 #
4 # Copyright (c) 2003-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) 2003-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 from esys.pycad import * #domain constructor
27 from esys.pycad.gmsh import Design #Finite Element meshing package
28 from esys.finley import MakeDomain #Converter for escript
29 import os #file path tool
30 import numpy as np #numerial python for arrays
31 from math import * # math package
32 #used to construct polygons for plotting from pycad
33 from cblib import getLoopCoords
34
35 #set modal to 1 for a syncline or -1 for an anticline structural
36 #configuration
37 modal=1
38
39 # the folder to put our outputs in, leave blank "" for script path -
40 # note this folder path must exist to work
41 save_path = "data/heatrefrac001"
42
43 # Overall Domain
44 p0=Point(0.0, 0.0, 0.0)
45 p1=Point(0.0, -6000.0, 0.0)
46 p2=Point(5000.0, -6000.0, 0.0)
47 p3=Point(5000.0, 0.0, 0.0)
48
49 l01=Line(p0, p1)
50 l12=Line(p1, p2)
51 l23=Line(p2, p3)
52 l30=Line(p3, p0)
53
54 c=CurveLoop(l01, l12, l23, l30)
55
56 # Material Boundary
57
58 x=[ Point(i*100.0,-2500+modal*1500.*cos(pi*i*100.0/2500.0+pi)) for i in range(0,51) ]
59 mysp = Spline(*tuple(x))
60 x1=Spline.getStartPoint(mysp)
61 x2=Spline.getEndPoint(mysp)
62
63 # TOP BLOCK
64 tbl1=Line(p0,x1)
65 tbl2=mysp
66 tbl3=Line(x2,p3)
67 tblockloop = CurveLoop(tbl1,tbl2,tbl3,l30)
68 tblock = PlaneSurface(tblockloop)
69
70 tpg = getLoopCoords(tblockloop)
71 np.savetxt(os.path.join(save_path,"toppg"),tpg,delimiter=" ")
72
73 # BOTTOM BLOCK
74 bbl1=Line(x1,p1)
75 bbl3=Line(p2,x2)
76 bbl4=-mysp
77 bblockloop = CurveLoop(bbl1,l12,bbl3,bbl4)
78 bblock = PlaneSurface(bblockloop)
79
80 #clockwise check as splines must be set as polygons in the point order
81 #they were created. Otherwise get a line across plot.
82 bblockloop2=CurveLoop(mysp,Line(x2,p2),Line(p2,p1),Line(p1,x1))
83 bpg = getLoopCoords(bblockloop2)
84 np.savetxt(os.path.join(save_path,"botpg"),bpg,delimiter=" ")
85
86 # Create a Design which can make the mesh
87 d=Design(dim=2, element_size=200)
88 # Add the subdomains and flux boundaries.
89 d.addItems(PropertySet("top",tblock),PropertySet("bottom",bblock),PropertySet("linebottom",l12))
90 # Create the geometry, mesh and Escript domain
91 d.setScriptFileName(os.path.join(save_path,"heatrefraction_mesh001.geo"))
92
93 d.setMeshFileName(os.path.join(save_path,"heatrefraction_mesh001.msh"))
94 domain=MakeDomain(d, integrationOrder=-1, reducedIntegrationOrder=-1, optimizeLabeling=True)
95 # Create a file that can be read back in to python with mesh=ReadMesh(fileName)
96 domain.write(os.path.join(save_path,"heatrefraction_mesh001.fly"))
97
98

  ViewVC Help
Powered by ViewVC 1.1.26