/[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 2648 - (show annotations)
Mon Sep 7 00:06:15 2009 UTC (11 years, 2 months ago) by jfenwick
File MIME type: text/x-python
File size: 3503 byte(s)
ecording some changes related to Anthony's cookbook.
These tests still do not run.

test_heatref.py is the current problem.
Current problems include:
undefined symbols
a dependence on X somehow.

I've already addressed the matplotlib X thing somewhere else.
Find and apply.


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 from cblib import needdirs
33
34
35 # routine to find consecutive coordinates of a loop in pycad
36 def loopcoords(loop):
37 # return all construction points of input
38 temp = loop.getConstructionPoints()
39 #create a numpy array for xyz components or construction points
40 coords = np.zeros([len(temp),3],float)
41 #place construction points in array
42 for i in range(0,len(temp)):
43 coords[i,:]=temp[i].getCoordinates()
44 #return a numpy array
45 return coords
46
47 save_path = "data/heatrefrac002"
48 needdirs([save_path])
49
50 # Overall Domain
51 p0=Point(0.0, 0.0, 0.0)
52 p1=Point(0.0, -6000.0, 0.0)
53 p2=Point(5000.0, -6000.0, 0.0)
54 p3=Point(5000.0, 0.0, 0.0)
55
56 l01=Line(p0, p1)
57 l12=Line(p1, p2)
58 l23=Line(p2, p3)
59 l30=Line(p3, p0)
60
61 c=CurveLoop(l01, l12, l23, l30)
62
63 # Material Boundary
64
65 p4=Point(0.0, -2400.0, 0.0)
66 p5=Point(2000.0, -2400.0, 0.0)
67 p6=Point(3000.0, -6000.0, 0.0)
68 p7=Point(5000.0, -2400.0, 0.0)
69
70 # TOP BLOCK
71 tbl1=Line(p0,p4)
72 tbl2=Line(p4,p5)
73 tbl3=Line(p5,p7)
74 tbl4=Line(p7,p3)
75 tblockloop = CurveLoop(tbl1,tbl2,tbl3,tbl4,l30)
76 tblock = PlaneSurface(tblockloop)
77
78 tpg = loopcoords(tblockloop)
79 np.savetxt(os.path.join(save_path,"toppg"),tpg,delimiter=" ")
80
81 # BOTTOM BLOCK LEFT
82 bbll1=Line(p4,p1)
83 bbll2=Line(p1,p6)
84 bbll3=Line(p6,p5)
85 bbll4=-tbl2
86 bblockloopl = CurveLoop(bbll1,bbll2,bbll3,bbll4)
87 bblockl = PlaneSurface(bblockloopl)
88
89 #clockwise check
90 #bblockloopl2=CurveLoop(mysp,Line(x2,p2),Line(p2,p1),Line(p1,x1))
91 bpg = loopcoords(bblockloopl)
92 np.savetxt(os.path.join(save_path,"botpgl"),bpg,delimiter=" ")
93
94 # BOTTOM BLOCK RIGHT
95 bbrl1=Line(p6,p2)
96 bbrl2=Line(p2,p7)
97 bbrl3=-tbl3
98 bbrl4=-bbll3
99 bblockloopr = CurveLoop(bbrl1,bbrl2,bbrl3,bbrl4)
100 bblockr = PlaneSurface(bblockloopr)
101
102 #clockwise check
103 #bblockloopr2=CurveLoop(mysp,Line(x2,p2),Line(p2,p1),Line(p1,x1))
104 bpg = loopcoords(bblockloopr)
105 np.savetxt(os.path.join(save_path,"botpgr"),bpg,delimiter=" ")
106
107 # Create a Design which can make the mesh
108 d=Design(dim=2, element_size=200)
109 # Add the subdomains and flux boundaries.
110 d.addItems(PropertySet("top",tblock),PropertySet("bottomleft",bblockl),PropertySet("bottomright",bblockr),PropertySet("linebottom",bbll21, bbrl1))
111 # Create the geometry, mesh and Escript domain
112 d.setScriptFileName(os.path.join(save_path,"heatrefraction_mesh003.geo"))
113
114 d.setMeshFileName(os.path.join(save_path,"heatrefraction_mesh003.msh"))
115 domain=MakeDomain(d, integrationOrder=-1, reducedIntegrationOrder=-1, optimizeLabeling=True)
116 # Create a file that can be read back in to python with mesh=ReadMesh(fileName)
117 domain.write(os.path.join(save_path,"heatrefraction_mesh003.fly"))
118

  ViewVC Help
Powered by ViewVC 1.1.26