/[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 3057 - (hide annotations)
Tue Jul 6 02:10:59 2010 UTC (11 years, 3 months ago) by ahallam
File MIME type: text/x-python
File size: 7564 byte(s)


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     # 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 ahallam 3057
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 ahallam 3055 #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 ahallam 3057 if (i > 0):
139 ahallam 3055 #print i,xs[i],ys[0],zs[i,0]
140 ahallam 3057 tp2=point_ar[0]; tp3=point_ar[-1]
141     tl0=Line(cp0,tp2); tl1=Line(cp1,tp3)
142 ahallam 3055 linex0_ar.append(tl0); linexy_ar.append(tl1)
143 ahallam 3057 cp0=point_ar[0]; cp1=point_ar[-1]
144 ahallam 3055
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 ahallam 3057 #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 ahallam 3055
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 ahallam 3057 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 ahallam 3055
184 ahallam 3057 #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 ahallam 3055
189 ahallam 3057 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 ahallam 3055 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 ahallam 3057 vintfa=Volume(SurfaceLoop(stop,*tuple(sintfa_ar+sintf_ar)))
200 ahallam 3055 vintfb=Volume(SurfaceLoop(sbot,*tuple(sintfb_ar+sintf_ar)))
201    
202 ahallam 3057 # Create the volume.
203     #sloop=SurfaceLoop(stop,sbot,*tuple(sintfa_ar+sintfb_ar))
204     #model=Volume(sloop)
205 ahallam 3055
206 ahallam 3057
207 ahallam 3055 #############################################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 ahallam 3057 #d.addItems(stop,sbot)
212 ahallam 3055 #d.addItems(PropertySet('intf',sintf))
213    
214 ahallam 3057 #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 ahallam 3055
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