/[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 5707 - (show annotations)
Mon Jun 29 03:59:06 2015 UTC (3 years, 8 months ago) by sshaw
File MIME type: text/x-python
File size: 6312 byte(s)
adding copyright headers to files without copyright info, moved header to top of file in some cases where it wasn't
1 ##############################################################################
2 #
3 # Copyright (c) 2009-2015 by The University of Queensland
4 # http://www.uq.edu.au
5 #
6 # Primary Business: Queensland, Australia
7 # Licensed under the Open Software License version 3.0
8 # http://www.opensource.org/licenses/osl-3.0.php
9 #
10 # Development until 2012 by Earth Systems Science Computational Center (ESSCC)
11 # Development 2012-2013 by School of Earth Sciences
12 # Development from 2014 by Centre for Geoscience Computing (GeoComp)
13 #
14 ##############################################################################
15 from __future__ import division, print_function
16
17 __copyright__="""Copyright (c) 2009-2015 by The University of Queensland
18 http://www.uq.edu.au
19 Primary Business: Queensland, Australia"""
20 __license__="""Licensed under the Open Software License version 3.0
21 http://www.opensource.org/licenses/osl-3.0.php"""
22 __url__="https://launchpad.net/escript-finley"
23
24 """
25 Author: Antony Hallam antony.hallam@uqconnect.edu.au
26 """
27
28
29 ############################################################FILE HEADER
30 # example09m.py
31 # Create a simple 3D model for use in example09. This is the low res
32 # mesh for illustration purposes only.
33 #
34 #######################################################EXTERNAL MODULES
35 from esys.pycad import * #domain constructor
36 from esys.pycad.extras import layer_cake
37 from esys.pycad.gmsh import Design #Finite Element meshing package
38 from esys.escript import mkDir, getMPISizeWorld
39 import os
40 import numpy as np
41 try:
42 # This imports the rectangle domain function
43 from esys.finley import MakeDomain
44 HAVE_FINLEY = True
45 except ImportError:
46 print("Finley module not available")
47 HAVE_FINLEY = False
48 ########################################################MPI WORLD CHECK
49 if getMPISizeWorld() > 1:
50 import sys
51 print("This example will not run in an MPI world.")
52 sys.exit(0)
53
54 if HAVE_FINLEY:
55 # make sure path exists
56 save_path= os.path.join("data","example09m")
57 mkDir(save_path)
58
59 ################################################ESTABLISHING PARAMETERS
60 # Time related variables.
61 testing=True
62 if testing:
63 print('This script is currently optioned for testing..')
64 print("Try changing the testing variable to False for more iterations.")
65 xwidth=40.
66 ywidth=40.
67 depth=20.
68 else:
69 #Model Parameters
70 xwidth=100.0 #x width of model
71 ywidth=100.0 #y width of model
72 depth=50.0 #depth of model
73
74 intf=depth/2. #Depth of the interface.
75 element_size=4.0
76
77 ####################################################DOMAIN CONSTRUCTION
78 # Domain Corners
79 p0=Point(0.0, 0.0, 0.0)
80 p1=Point(xwidth, 0.0, 0.0)
81 p2=Point(xwidth, ywidth, 0.0)
82 p3=Point(0.0, ywidth, 0.0)
83 p4=Point(0.0, ywidth, depth)
84 p5=Point(0.0, 0.0, depth)
85 p6=Point(xwidth, 0.0, depth)
86 p7=Point(xwidth, ywidth, depth)
87 # Join corners in anti-clockwise manner.
88 l01=Line(p0, p1)
89 l12=Line(p1, p2)
90 l23=Line(p2, p3)
91 l30=Line(p3, p0)
92 l56=Line(p5, p6)
93 l67=Line(p6, p7)
94 l74=Line(p7, p4)
95 l45=Line(p4, p5)
96
97 # Join line segments to create domain boundaries and then surfaces
98 ctop=CurveLoop(l01, l12, l23, l30); stop=PlaneSurface(ctop)
99 cbot=CurveLoop(-l67, -l56, -l45, -l74); sbot=PlaneSurface(cbot)
100
101 # for each side
102 ip0=Point(0.0, 0.0, intf)
103 ip1=Point(xwidth, 0.0, intf)
104 ip2=Point(xwidth, ywidth, intf)
105 ip3=Point(0.0, ywidth, intf)
106
107 linte_ar=[]; #lines for vertical edges
108 linhe_ar=[]; #lines for horizontal edges
109 linte_ar.append(Line(p0,ip0))
110 linte_ar.append(Line(ip0,p5))
111 linte_ar.append(Line(p1,ip1))
112 linte_ar.append(Line(ip1,p6))
113 linte_ar.append(Line(p2,ip2))
114 linte_ar.append(Line(ip2,p7))
115 linte_ar.append(Line(p3,ip3))
116 linte_ar.append(Line(ip3,p4))
117
118 linhe_ar.append(Line(ip0,ip1))
119 linhe_ar.append(Line(ip1,ip2))
120 linhe_ar.append(Line(ip2,ip3))
121 linhe_ar.append(Line(ip3,ip0))
122
123 cintfa_ar=[]; cintfb_ar=[] #curveloops for above and below interface on sides
124 cintfa_ar.append(CurveLoop(linte_ar[0],linhe_ar[0],-linte_ar[2],-l01))
125 cintfa_ar.append(CurveLoop(linte_ar[2],linhe_ar[1],-linte_ar[4],-l12))
126 cintfa_ar.append(CurveLoop(linte_ar[4],linhe_ar[2],-linte_ar[6],-l23))
127 cintfa_ar.append(CurveLoop(linte_ar[6],linhe_ar[3],-linte_ar[0],-l30))
128
129 cintfb_ar.append(CurveLoop(linte_ar[1],l56,-linte_ar[3],-linhe_ar[0]))
130 cintfb_ar.append(CurveLoop(linte_ar[3],l67,-linte_ar[5],-linhe_ar[1]))
131 cintfb_ar.append(CurveLoop(linte_ar[5],l74,-linte_ar[7],-linhe_ar[2]))
132 cintfb_ar.append(CurveLoop(linte_ar[7],l45,-linte_ar[1],-linhe_ar[3]))
133
134 sintfa_ar=[PlaneSurface(cintfa_ar[i]) for i in range(0,4)]
135 sintfb_ar=[PlaneSurface(cintfb_ar[i]) for i in range(0,4)]
136
137 sintf=PlaneSurface(CurveLoop(*tuple(linhe_ar)))
138
139 vintfa=Volume(SurfaceLoop(stop,-sintf,*tuple(sintfa_ar)))
140 vintfb=Volume(SurfaceLoop(sbot,sintf,*tuple(sintfb_ar)))
141
142 #############################################EXPORTING MESH FOR ESCRIPT
143 # Create a Design which can make the mesh
144 d=Design(dim=3, element_size=element_size, order=2)
145
146 d.addItems(PropertySet('vintfa',vintfa))
147 d.addItems(PropertySet('vintfb',vintfb))
148 d.addItems(PropertySet('stop',stop))
149 d.addItems(PropertySet('sbot',sbot))
150
151 d.setScriptFileName(os.path.join(save_path,"example09m.geo"))
152 d.setMeshFileName(os.path.join(save_path,"example09m.msh"))
153 #
154 # make the domain:
155 #
156 domain=MakeDomain(d)
157 # Create a file that can be read back in to python with
158 # mesh=ReadMesh(fileName)
159 domain.write(os.path.join(save_path,"example09m.fly"))
160
161 if testing:
162 intfaces=np.array([10,30,50,55,80,100,200,250,400])/100.
163 else:
164 intfaces=np.array([10,30,50,55,80,100,200,250,400])/10.
165
166 # Specify the domain.
167 domaindes=Design(dim=3,element_size=element_size,order=2)
168 cmplx_domain=layer_cake(domaindes,xwidth,ywidth,intfaces)
169 cmplx_domain.setScriptFileName(os.path.join(save_path,"example09lc.geo"))
170 cmplx_domain.setMeshFileName(os.path.join(save_path,"example09lc.msh"))
171 dcmplx=MakeDomain(cmplx_domain)
172 dcmplx.write(os.path.join(save_path,"example09lc.fly"))
173
174
175
176

  ViewVC Help
Powered by ViewVC 1.1.26