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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 2597 - (show annotations)
Thu Aug 6 03:30:09 2009 UTC (10 years, 3 months ago) by ahallam
File MIME type: text/x-python
File size: 3453 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 *
27 from esys.pycad.gmsh import Design
28 from esys.finley import MakeDomain
29 import os
30 import numpy as np
31 from math import *
32
33 # routine to find consecutive coordinates of a loop in pycad
34 def loopcoords(loop):
35 # return all construction points of input
36 temp = loop.getConstructionPoints()
37 #create a numpy array for xyz components or construction points
38 coords = np.zeros([len(temp),3],float)
39 #place construction points in array
40 for i in range(0,len(temp)):
41 coords[i,:]=temp[i].getCoordinates()
42 #return a numpy array
43 return coords
44
45 save_path = "data/heatrefrac002"
46
47 # Overall Domain
48 p0=Point(0.0, 0.0, 0.0)
49 p1=Point(0.0, -6000.0, 0.0)
50 p2=Point(5000.0, -6000.0, 0.0)
51 p3=Point(5000.0, 0.0, 0.0)
52
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 # Material Boundary
61
62 p4=Point(0.0, -2400.0, 0.0)
63 p5=Point(2000.0, -2400.0, 0.0)
64 p6=Point(3000.0, -6000.0, 0.0)
65 p7=Point(5000.0, -2400.0, 0.0)
66
67 # TOP BLOCK
68 tbl1=Line(p0,p4)
69 tbl2=Line(p4,p5)
70 tbl3=Line(p5,p7)
71 tbl4=Line(p7,p3)
72 tblockloop = CurveLoop(tbl1,tbl2,tbl3,tbl4,l30)
73 tblock = PlaneSurface(tblockloop)
74
75 tpg = loopcoords(tblockloop)
76 np.savetxt(os.path.join(save_path,"toppg"),tpg,delimiter=" ")
77
78 # BOTTOM BLOCK LEFT
79 bbll1=Line(p4,p1)
80 bbll2=Line(p1,p6)
81 bbll3=Line(p6,p5)
82 bbll4=-tbl2
83 bblockloopl = CurveLoop(bbll1,bbll2,bbll3,bbll4)
84 bblockl = PlaneSurface(bblockloopl)
85
86 #clockwise check
87 #bblockloopl2=CurveLoop(mysp,Line(x2,p2),Line(p2,p1),Line(p1,x1))
88 bpg = loopcoords(bblockloopl)
89 np.savetxt(os.path.join(save_path,"botpgl"),bpg,delimiter=" ")
90
91 # BOTTOM BLOCK RIGHT
92 bbrl1=Line(p6,p2)
93 bbrl2=Line(p2,p7)
94 bbrl3=-tbl3
95 bbrl4=-bbll3
96 bblockloopr = CurveLoop(bbrl1,bbrl2,bbrl3,bbrl4)
97 bblockr = PlaneSurface(bblockloopr)
98
99 #clockwise check
100 #bblockloopr2=CurveLoop(mysp,Line(x2,p2),Line(p2,p1),Line(p1,x1))
101 bpg = loopcoords(bblockloopr)
102 np.savetxt(os.path.join(save_path,"botpgr"),bpg,delimiter=" ")
103
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),PropertySet("bottomleft",bblockl),PropertySet("bottomright",bblockr),PropertySet("linebottom",bbll21, bbrl1))
108 # Create the geometry, mesh and Escript domain
109 d.setScriptFileName(os.path.join(save_path,"heatrefraction_mesh003.geo"))
110
111 d.setMeshFileName(os.path.join(save_path,"heatrefraction_mesh003.msh"))
112 domain=MakeDomain(d, integrationOrder=-1, reducedIntegrationOrder=-1, optimizeLabeling=True)
113 # Create a file that can be read back in to python with mesh=ReadMesh(fileName)
114 domain.write(os.path.join(save_path,"heatrefraction_mesh003.fly"))
115

  ViewVC Help
Powered by ViewVC 1.1.26