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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 3057 - (show annotations)
Tue Jul 6 02:10:59 2010 UTC (11 years, 5 months ago) by ahallam
File MIME type: text/x-python
File size: 7564 byte(s)


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
26
27 ############################################################FILE HEADER
28 # example09m.py
29 # Create a simple 3D model for use in example09.
30 #
31 #######################################################EXTERNAL MODULES
32 from esys.pycad import * #domain constructor
33 from esys.pycad.gmsh import Design #Finite Element meshing package
34 from esys.finley import MakeDomain #Converter for escript
35 from esys.escript import mkDir, getMPISizeWorld
36 from esys.escript.unitsSI import *
37 import os
38 from math import *
39 import pylab as pl
40 import numpy as np
41 ########################################################MPI WORLD CHECK
42 if getMPISizeWorld() > 1:
43 import sys
44 print "This example will not run in an MPI world."
45 sys.exit(0)
46
47 # make sure path exists
48 save_path= os.path.join("data","example09")
49 mkDir(save_path)
50
51 ################################################ESTABLISHING PARAMETERS
52 #Model Parameters
53 xwidth=2000.0*m #x width of model
54 ywidth=2000.0*m #y width of model
55 depth=500.0*m #depth of model
56
57 ####################################################DOMAIN CONSTRUCTION
58 # Domain Corners
59 p0=Point(0.0, 0.0, 0.0)
60 p1=Point(xwidth, 0.0, 0.0)
61 p2=Point(xwidth, ywidth, 0.0)
62 p3=Point(0.0, ywidth, 0.0)
63 p4=Point(0.0, ywidth, depth)
64 p5=Point(0.0, 0.0, depth)
65 p6=Point(xwidth, 0.0, depth)
66 p7=Point(xwidth, ywidth, depth)
67 # Join corners in anti-clockwise manner.
68 l01=Line(p0, p1)
69 l12=Line(p1, p2)
70 l23=Line(p2, p3)
71 l30=Line(p3, p0)
72 l56=Line(p5, p6)
73 l67=Line(p6, p7)
74 l74=Line(p7, p4)
75 l45=Line(p4, p5)
76
77 # Join line segments to create domain boundaries and then surfaces
78 ctop=CurveLoop(l01, l12, l23, l30); stop=PlaneSurface(ctop)
79 cbot=CurveLoop(-l67, -l56, -l45, -l74); sbot=PlaneSurface(cbot)
80
81 #create interface
82 sspl=51 #number of discrete points in spline
83 xdsp=xwidth/(sspl-1) #dx of spline steps for xwidth
84 ydsp=ywidth/(sspl-1) #dy of spline steps for ywidth
85 dep_sp=250.0*m #avg depth of spline
86 h_sp=25.0*m #heigh of spline
87 orit=-1.0 #orientation of spline 1.0=>up -1.0=>down
88
89 # Generate Material Boundary
90 #x=[ Point(i*xdsp\
91 # ,-dep_sp+orit*h_sp*cos(pi*i*xdsp/dep_sp+pi))\
92 # for i in range(0,sspl)\
93 #]
94 #function along x
95 x=[1+orit*cos(2*sspl*pi*i*xdsp/xwidth)\
96 for i in range(0,sspl)]
97 #function along y
98 y=[1+orit*cos(2*sspl*pi*i*ydsp/ywidth)\
99 for i in range(0,sspl)]
100 #containers for surface data
101 xs=np.zeros(sspl,'float')
102 ys=np.zeros(sspl,'float')
103 zs=np.zeros([sspl,sspl],'float')
104
105 for i in range(0,sspl):
106 xs[i]=i*xdsp
107 for j in range(0,sspl):
108 if (i == 0):
109 ys[j]=j*ydsp
110 zs[i,j]=x[i]*y[j]*h_sp+dep_sp
111
112 pl.plot(x); pl.plot(y,'ro')
113 pl.savefig(os.path.join(save_path,"interface.png"))
114 pl.clf()
115 pl.contourf(xs,ys,zs)
116 pl.savefig(os.path.join(save_path,"interface_cf.png"))
117
118 spl_ar=[] #interface spline array
119 linex0_ar=[] #interface line array along x
120 linexy_ar=[] #interface line array along x,ymax
121 sintf_ar=[] #interface surface array
122 nsintf_ar=[]
123 #loop through x
124 for i in range(0,sspl):
125 #create all the points required
126 point_ar=[Point(xs[i],ys[j],zs[i,j]) for j in range(0,sspl)]
127
128 if (i == 0):
129 ip0=point_ar[0]
130 ip3=point_ar[-1]
131 if (i == sspl-1):
132 ip1=point_ar[0]
133 ip2=point_ar[-1]
134
135 #create the spline and store it
136 spl_ar.append(Spline(*tuple(point_ar)))
137 #create the lines along the x axis and x axis at ymax
138 if (i > 0):
139 #print i,xs[i],ys[0],zs[i,0]
140 tp2=point_ar[0]; tp3=point_ar[-1]
141 tl0=Line(cp0,tp2); tl1=Line(cp1,tp3)
142 linex0_ar.append(tl0); linexy_ar.append(tl1)
143 cp0=point_ar[0]; cp1=point_ar[-1]
144
145 for i in range(0,sspl):
146 #create the mini surfaces via curveloops and then ruledsurfaces
147 if (i < sspl-1):
148 #print 'i',i
149 tc0=CurveLoop(linex0_ar[i],spl_ar[i+1],-linexy_ar[i],-spl_ar[i])
150 ntc0=CurveLoop(spl_ar[i],linexy_ar[i],-spl_ar[i+1],-linex0_ar[i])
151 ts0=RuledSurface(tc0); sintf_ar.append(ts0)
152 nts0=RuledSurface(ntc0); nsintf_ar.append(nts0)
153 #create the interface using the surface loop constructor
154 sintf=SurfaceLoop(*tuple(sintf_ar))
155
156 # for each side
157 #ip0=Point(xs[0],ys[0],zs[0,0])
158 #ip1=Point(xs[-1],ys[0],zs[-1,0])
159 #ip2=Point(xs[-1],ys[-1],zs[-1,-1])
160 #ip3=Point(xs[0],ys[-1],zs[0,-1])
161
162 linte_ar=[]; #lines for vertical edges
163 linhe_ar=[]; #lines for horizontal edges
164 linte_ar.append(Line(p0,ip0))
165 linte_ar.append(Line(ip0,p5))
166 linte_ar.append(Line(p1,ip1))
167 linte_ar.append(Line(ip1,p6))
168 linte_ar.append(Line(p2,ip2))
169 linte_ar.append(Line(ip2,p7))
170 linte_ar.append(Line(p3,ip3))
171 linte_ar.append(Line(ip3,p4))
172
173 linhe_ar.append(Line(ip0,ip1))
174 linhe_ar.append(Line(ip1,ip2))
175 linhe_ar.append(Line(ip2,ip3))
176 linhe_ar.append(Line(ip3,ip0))
177
178 cintfa_ar=[]; cintfb_ar=[] #curveloops for above and below interface on sides
179 cintfa_ar.append(CurveLoop(-linte_ar[2],-l01,linte_ar[0],*tuple(linex0_ar)))
180 cintfa_ar.append(CurveLoop(linte_ar[2],spl_ar[-1],-linte_ar[4],-l12))
181 cintfa_ar.append(CurveLoop(-linte_ar[6],-l23,linte_ar[4],*tuple(linexy_ar)))
182 cintfa_ar.append(CurveLoop(linte_ar[6],-spl_ar[0],-linte_ar[0],-l30))
183
184 #cintfb_ar.append(CurveLoop(linte_ar[1],l56,-linte_ar[3],-linhe_ar[0]))
185 #cintfb_ar.append(CurveLoop(linte_ar[3],l67,-linte_ar[5],-linhe_ar[1]))
186 #cintfb_ar.append(CurveLoop(linte_ar[5],l74,-linte_ar[7],-linhe_ar[2]))
187 #cintfb_ar.append(CurveLoop(linte_ar[7],l45,-linte_ar[1],-linhe_ar[3]))
188
189 cintfb_ar.append(-CurveLoop(linte_ar[3],-l56,-linte_ar[1],*tuple(linex0_ar)))
190 cintfb_ar.append(-CurveLoop(linte_ar[5],-l67,-linte_ar[3],-spl_ar[-1]))
191 cintfb_ar.append(CurveLoop(linte_ar[5],l74,-linte_ar[7],*tuple(linexy_ar)))
192 cintfb_ar.append(CurveLoop(linte_ar[7],l45,-linte_ar[1],spl_ar[0]))
193
194
195
196 sintfa_ar=[PlaneSurface(cintfa_ar[i]) for i in range(0,4)]
197 sintfb_ar=[PlaneSurface(cintfb_ar[i]) for i in range(0,4)]
198
199 vintfa=Volume(SurfaceLoop(stop,*tuple(sintfa_ar+sintf_ar)))
200 vintfb=Volume(SurfaceLoop(sbot,*tuple(sintfb_ar+sintf_ar)))
201
202 # Create the volume.
203 #sloop=SurfaceLoop(stop,sbot,*tuple(sintfa_ar+sintfb_ar))
204 #model=Volume(sloop)
205
206
207 #############################################EXPORTING MESH FOR ESCRIPT
208 # Create a Design which can make the mesh
209 d=Design(dim=3, element_size=200.0*m)
210 # Add the subdomains and flux boundaries.
211 #d.addItems(stop,sbot)
212 #d.addItems(PropertySet('intf',sintf))
213
214 #d.addItems(PropertySet('vintfa',vintfa))
215 #d.addItems(PropertySet('vintfb',vintfb))
216 #d.addItems(PropertySet('l1',linte_ar[0]))
217 d.addItems(*tuple(sintfb_ar))
218 d.addItems(sbot)
219 d.addItems(*tuple(sintf_ar))
220
221 d.setScriptFileName(os.path.join(save_path,"example09m.geo"))
222
223 #d.setMeshFileName(os.path.join(save_path,"example09m.msh"))
224 #
225 # make the finley domain:
226 #
227 domain=MakeDomain(d)
228 # Create a file that can be read back in to python with
229 # mesh=ReadMesh(fileName)
230 #domain.write(os.path.join(save_path,"example09m.fly"))
231
232

  ViewVC Help
Powered by ViewVC 1.1.26