/[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 3055 - (show annotations)
Sun Jul 4 21:52:36 2010 UTC (11 years, 3 months ago) by ahallam
File MIME type: text/x-python
File size: 7396 byte(s)
3d wave equation and ABC
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 l05=Line(p0, p5)
77 l16=Line(p1, p6)
78 l27=Line(p2, p7)
79 l34=Line(p3, p4)
80
81 # Join line segments to create domain boundaries and then surfaces
82 ctop=CurveLoop(l01, l12, l23, l30); stop=PlaneSurface(ctop)
83 cbot=CurveLoop(-l67, -l56, -l45, -l74); sbot=PlaneSurface(cbot)
84 cx0 =CurveLoop(l05, l56, -l16, -l01); sx0=PlaneSurface(cx0)
85 cy0 =CurveLoop(-l30, l34, l45, -l05); sy0=PlaneSurface(cy0)
86 cxy =CurveLoop(-l34, -l23, l27, l74); sxy=PlaneSurface(cxy)
87 cyx =CurveLoop(-l12, l16, l67, -l27); syx=PlaneSurface(cyx)
88
89 # Create the volume.
90 sloop=SurfaceLoop(stop,sbot,sx0,sy0,sxy,syx)
91 model=Volume(sloop)
92
93
94 #create interface
95 sspl=51 #number of discrete points in spline
96 xdsp=xwidth/(sspl-1) #dx of spline steps for xwidth
97 ydsp=ywidth/(sspl-1) #dy of spline steps for ywidth
98 dep_sp=250.0*m #avg depth of spline
99 h_sp=25.0*m #heigh of spline
100 orit=-1.0 #orientation of spline 1.0=>up -1.0=>down
101
102 # Generate Material Boundary
103 #x=[ Point(i*xdsp\
104 # ,-dep_sp+orit*h_sp*cos(pi*i*xdsp/dep_sp+pi))\
105 # for i in range(0,sspl)\
106 #]
107 #function along x
108 x=[1+orit*cos(2*sspl*pi*i*xdsp/xwidth)\
109 for i in range(0,sspl)]
110 #function along y
111 y=[1+orit*cos(2*sspl*pi*i*ydsp/ywidth)\
112 for i in range(0,sspl)]
113 #containers for surface data
114 xs=np.zeros(sspl,'float')
115 ys=np.zeros(sspl,'float')
116 zs=np.zeros([sspl,sspl],'float')
117
118 for i in range(0,sspl):
119 xs[i]=i*xdsp
120 for j in range(0,sspl):
121 if (i == 0):
122 ys[j]=j*ydsp
123 zs[i,j]=x[i]*y[j]*h_sp+dep_sp
124
125 pl.plot(x); pl.plot(y,'ro')
126 pl.savefig(os.path.join(save_path,"interface.png"))
127 pl.clf()
128 pl.contourf(xs,ys,zs)
129 pl.savefig(os.path.join(save_path,"interface_cf.png"))
130
131 spl_ar=[] #interface spline array
132 linex0_ar=[] #interface line array along x
133 linexy_ar=[] #interface line array along x,ymax
134 sintf_ar=[] #interface surface array
135 nsintf_ar=[]
136 #loop through x
137 for i in range(0,sspl):
138 #create all the points required
139 point_ar=[Point(xs[i],ys[j],zs[i,j]) for j in range(0,sspl)]
140 #create the spline and store it
141 spl_ar.append(Spline(*tuple(point_ar)))
142 #create the lines along the x axis and x axis at ymax
143 if (i < sspl-1):
144 #print i,xs[i],ys[0],zs[i,0]
145 tp0=Point(xs[i],ys[0],zs[i,0])
146 tp1=Point(xs[i+1],ys[0],zs[i+1,0])
147 tp2=Point(xs[i],ys[-1],zs[i,-1])
148 tp3=Point(xs[i+1],ys[-1],zs[i+1,-1])
149 tl0=Line(tp0,tp1); tl1=Line(tp2,tp3)
150 linex0_ar.append(tl0); linexy_ar.append(tl1)
151
152 for i in range(0,sspl):
153 #create the mini surfaces via curveloops and then ruledsurfaces
154 if (i < sspl-1):
155 #print 'i',i
156 tc0=CurveLoop(linex0_ar[i],spl_ar[i+1],-linexy_ar[i],-spl_ar[i])
157 ntc0=CurveLoop(spl_ar[i],linexy_ar[i],-spl_ar[i+1],-linex0_ar[i])
158 ts0=RuledSurface(tc0); sintf_ar.append(ts0)
159 nts0=RuledSurface(ntc0); nsintf_ar.append(nts0)
160 #create the interface using the surface loop constructor
161 sintf=SurfaceLoop(*tuple(sintf_ar))
162
163 # for each side
164 ip0=Point(xs[0],ys[0],zs[0,0])
165 ip1=Point(xs[-1],ys[0],zs[-1,0])
166 ip2=Point(xs[-1],ys[-1],zs[-1,-1])
167 ip3=Point(xs[0],ys[-1],zs[0,-1])
168
169 linte_ar=[]; #lines for vertical edges
170 linhe_ar=[]; #lines for horizontal edges
171 linte_ar.append(Line(p0,ip0))
172 linte_ar.append(Line(ip0,p5))
173 linte_ar.append(Line(p1,ip1))
174 linte_ar.append(Line(ip1,p6))
175 linte_ar.append(Line(p2,ip2))
176 linte_ar.append(Line(ip2,p7))
177 linte_ar.append(Line(p3,ip3))
178 linte_ar.append(Line(ip3,p4))
179
180 linhe_ar.append(Line(ip0,ip1))
181 linhe_ar.append(Line(ip1,ip2))
182 linhe_ar.append(Line(ip2,ip3))
183 linhe_ar.append(Line(ip3,ip0))
184
185 cintfa_ar=[]; cintfb_ar=[] #curveloops for above and below interface on sides
186 cintfa_ar.append(CurveLoop(linte_ar[0],linhe_ar[0],-linte_ar[2],-l01))
187 cintfa_ar.append(CurveLoop(linte_ar[2],linhe_ar[1],-linte_ar[4],-l12))
188 cintfa_ar.append(CurveLoop(linte_ar[4],linhe_ar[2],-linte_ar[6],-l23))
189 cintfa_ar.append(CurveLoop(linte_ar[6],linhe_ar[3],-linte_ar[0],-l30))
190
191 cintfb_ar.append(CurveLoop(linte_ar[1],l56,-linte_ar[3],-linhe_ar[0]))
192 cintfb_ar.append(CurveLoop(linte_ar[3],l67,-linte_ar[5],-linhe_ar[1]))
193 cintfb_ar.append(CurveLoop(linte_ar[5],l74,-linte_ar[7],-linhe_ar[2]))
194 cintfb_ar.append(CurveLoop(linte_ar[7],l45,-linte_ar[1],-linhe_ar[3]))
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+nsintf_ar)))
200 vintfb=Volume(SurfaceLoop(sbot,*tuple(sintfb_ar+sintf_ar)))
201
202
203 #############################################EXPORTING MESH FOR ESCRIPT
204 # Create a Design which can make the mesh
205 d=Design(dim=3, element_size=200.0*m)
206 # Add the subdomains and flux boundaries.
207 d.addItems(model,stop,sbot,sx0,sy0,sxy,syx)
208 #d.addItems(PropertySet('intf',sintf))
209
210 d.addItems(PropertySet('vintfa',vintfa),PropertySet('vintfb',vintfb))
211
212
213 d.setScriptFileName(os.path.join(save_path,"example09m.geo"))
214
215 #d.setMeshFileName(os.path.join(save_path,"example09m.msh"))
216 #
217 # make the finley domain:
218 #
219 domain=MakeDomain(d)
220 # Create a file that can be read back in to python with
221 # mesh=ReadMesh(fileName)
222 #domain.write(os.path.join(save_path,"example09m.fly"))
223
224

  ViewVC Help
Powered by ViewVC 1.1.26