/[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 3055 - (hide 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 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     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