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

Annotation of /trunk/doc/examples/cookbook/example09m.py

Parent Directory Parent Directory | Revision Log Revision Log


Revision 3069 - (hide annotations)
Wed Jul 21 03:24:48 2010 UTC (11 years, 3 months ago) by ahallam
File MIME type: text/x-python
File size: 5278 byte(s)
Updates to layer cake and cookbook as well as examples.
1 ahallam 3055
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    
26    
27     ############################################################FILE HEADER
28     # example09m.py
29 ahallam 3069 # Create a simple 3D model for use in example09. This is the low res
30     # mesh for illustration purposes only.
31 ahallam 3055 #
32     #######################################################EXTERNAL MODULES
33     from esys.pycad import * #domain constructor
34     from esys.pycad.gmsh import Design #Finite Element meshing package
35     from esys.finley import MakeDomain #Converter for escript
36     from esys.escript import mkDir, getMPISizeWorld
37     import os
38     ########################################################MPI WORLD CHECK
39     if getMPISizeWorld() > 1:
40     import sys
41     print "This example will not run in an MPI world."
42     sys.exit(0)
43    
44     # make sure path exists
45     save_path= os.path.join("data","example09")
46     mkDir(save_path)
47    
48     ################################################ESTABLISHING PARAMETERS
49     #Model Parameters
50 ahallam 3069 xwidth=1000.0 #x width of model
51     ywidth=1000.0 #y width of model
52     depth=200.0 #depth of model
53     intf=depth/2. #Depth of the interface.
54     element_size=5.0
55 ahallam 3055
56     ####################################################DOMAIN CONSTRUCTION
57     # Domain Corners
58     p0=Point(0.0, 0.0, 0.0)
59     p1=Point(xwidth, 0.0, 0.0)
60     p2=Point(xwidth, ywidth, 0.0)
61     p3=Point(0.0, ywidth, 0.0)
62     p4=Point(0.0, ywidth, depth)
63     p5=Point(0.0, 0.0, depth)
64     p6=Point(xwidth, 0.0, depth)
65     p7=Point(xwidth, ywidth, depth)
66     # Join corners in anti-clockwise manner.
67     l01=Line(p0, p1)
68     l12=Line(p1, p2)
69     l23=Line(p2, p3)
70     l30=Line(p3, p0)
71     l56=Line(p5, p6)
72     l67=Line(p6, p7)
73     l74=Line(p7, p4)
74     l45=Line(p4, p5)
75    
76     # Join line segments to create domain boundaries and then surfaces
77     ctop=CurveLoop(l01, l12, l23, l30); stop=PlaneSurface(ctop)
78     cbot=CurveLoop(-l67, -l56, -l45, -l74); sbot=PlaneSurface(cbot)
79    
80     # for each side
81 ahallam 3069 ip0=Point(0.0, 0.0, intf)
82     ip1=Point(xwidth, 0.0, intf)
83     ip2=Point(xwidth, ywidth, intf)
84     ip3=Point(0.0, ywidth, intf)
85 ahallam 3055
86     linte_ar=[]; #lines for vertical edges
87     linhe_ar=[]; #lines for horizontal edges
88     linte_ar.append(Line(p0,ip0))
89     linte_ar.append(Line(ip0,p5))
90     linte_ar.append(Line(p1,ip1))
91     linte_ar.append(Line(ip1,p6))
92     linte_ar.append(Line(p2,ip2))
93     linte_ar.append(Line(ip2,p7))
94     linte_ar.append(Line(p3,ip3))
95     linte_ar.append(Line(ip3,p4))
96    
97     linhe_ar.append(Line(ip0,ip1))
98     linhe_ar.append(Line(ip1,ip2))
99     linhe_ar.append(Line(ip2,ip3))
100     linhe_ar.append(Line(ip3,ip0))
101    
102     cintfa_ar=[]; cintfb_ar=[] #curveloops for above and below interface on sides
103 ahallam 3069 cintfa_ar.append(CurveLoop(linte_ar[0],linhe_ar[0],-linte_ar[2],-l01))
104     cintfa_ar.append(CurveLoop(linte_ar[2],linhe_ar[1],-linte_ar[4],-l12))
105     cintfa_ar.append(CurveLoop(linte_ar[4],linhe_ar[2],-linte_ar[6],-l23))
106     cintfa_ar.append(CurveLoop(linte_ar[6],linhe_ar[3],-linte_ar[0],-l30))
107 ahallam 3055
108 ahallam 3069 cintfb_ar.append(CurveLoop(linte_ar[1],l56,-linte_ar[3],-linhe_ar[0]))
109     cintfb_ar.append(CurveLoop(linte_ar[3],l67,-linte_ar[5],-linhe_ar[1]))
110     cintfb_ar.append(CurveLoop(linte_ar[5],l74,-linte_ar[7],-linhe_ar[2]))
111     cintfb_ar.append(CurveLoop(linte_ar[7],l45,-linte_ar[1],-linhe_ar[3]))
112 ahallam 3055
113     sintfa_ar=[PlaneSurface(cintfa_ar[i]) for i in range(0,4)]
114     sintfb_ar=[PlaneSurface(cintfb_ar[i]) for i in range(0,4)]
115    
116 ahallam 3069 sintf=PlaneSurface(CurveLoop(*tuple(linhe_ar)))
117 ahallam 3055
118 ahallam 3069 vintfa=Volume(SurfaceLoop(stop,-sintf,*tuple(sintfa_ar)))
119     vintfb=Volume(SurfaceLoop(sbot,sintf,*tuple(sintfb_ar)))
120    
121 ahallam 3057 # Create the volume.
122     #sloop=SurfaceLoop(stop,sbot,*tuple(sintfa_ar+sintfb_ar))
123     #model=Volume(sloop)
124 ahallam 3055
125 ahallam 3057
126 ahallam 3055 #############################################EXPORTING MESH FOR ESCRIPT
127     # Create a Design which can make the mesh
128 ahallam 3069 d=Design(dim=3, element_size=5.0, order=2)
129 ahallam 3055
130 ahallam 3069 d.addItems(PropertySet('vintfa',vintfa))
131     d.addItems(PropertySet('vintfb',vintfb))
132     d.addItems(PropertySet('stop',stop))
133     d.addItems(PropertySet('sbpt',sbot))
134 ahallam 3055
135     d.setScriptFileName(os.path.join(save_path,"example09m.geo"))
136 ahallam 3069 d.setMeshFileName(os.path.join(save_path,"example09m.msh"))
137 ahallam 3055 #
138     # make the finley domain:
139     #
140     domain=MakeDomain(d)
141     # Create a file that can be read back in to python with
142     # mesh=ReadMesh(fileName)
143 ahallam 3069 domain.write(os.path.join(save_path,"example09m.fly"))
144 ahallam 3055
145 ahallam 3069 from esys.pycad import layer_cake
146     intfaces=[10,30,50,55,80,100,200,250,400,500]
147 ahallam 3055
148 ahallam 3069 # Specify the domain.
149     domaindes=Design(dim=3,element_size=element_size,order=2)
150     cmplx_domain=layer_cake.LayerCake(domaindes,xwidth,ywidth,intfaces)
151     cmplx_domain.setScriptFileName(os.path.join(save_path,"example09lc.geo"))
152     cmplx_domain.setMeshFileName(os.path.join(save_path,"example09lc.msh"))
153     dcmplx=MakeDomain(cmplx_domain)
154     dcmplx.write(os.path.join(save_path,"example09lc.fly"))
155    
156    
157    
158    

  ViewVC Help
Powered by ViewVC 1.1.26