/[escript]/trunk/pycad/py_src/extras.py
ViewVC logotype

Contents of /trunk/pycad/py_src/extras.py

Parent Directory Parent Directory | Revision Log Revision Log


Revision 3062 - (show annotations)
Wed Jul 14 00:18:30 2010 UTC (9 years, 4 months ago) by ahallam
File MIME type: text/x-python
File size: 6048 byte(s)
Joels suggested changes for naming conventions. The function is now located in esys.pycad.extras as layer_cake.
1 # -*- coding: utf-8 -*-
2
3 ########################################################
4 #
5 # Copyright (c) 2003-2010 by University of Queensland
6 # Earth Systems Science Computational Center (ESSCC)
7 # http://www.uq.edu.au/esscc
8 #
9 # Primary Business: Queensland, Australia
10 # Licensed under the Open Software License version 3.0
11 # http://www.opensource.org/licenses/osl-3.0.php
12 #
13 ########################################################
14
15 __copyright__="""Copyright (c) 2003-2010 by University of Queensland
16 Earth Systems Science Computational Center (ESSCC)
17 http://www.uq.edu.au/esscc
18 Primary Business: Queensland, Australia"""
19 __license__="""Licensed under the Open Software License version 3.0
20 http://www.opensource.org/licenses/osl-3.0.php"""
21 __url__="https://launchpad.net/escript-finley"
22
23 """
24 Layer Cake
25
26 This script uses the pycad module to build a rectangular layered
27 system of prisms for modelling.
28
29 :var __author__: name of author
30 :var __copyright__: copyrights
31 :var __license__: licence agreement
32 :var __url__: url entry point on documentation
33 :var __version__: version
34 :var __date__: date of the version
35 """
36
37 __author__="Antony Hallam antony.hallam@uqconnect.edu.au"
38
39 from esys.pycad import * #domain constructor
40 from esys.pycad.gmsh import Design #Finite Element meshing package
41
42 def buildFreeSurface(xwidth,ywidth):
43 '''
44 Build a free surface to start the layer cake model.
45 This surface is planar.
46 Parameters:
47 xwidth :: width in x direction in meters.
48 ywidth :: width in y direction in meters.
49 '''
50
51 # Layer Corners
52 corner_points=[]
53 corner_points.append(Point(0.0, 0.0, 0.0))
54 corner_points.append(Point(xwidth, 0.0, 0.0))
55 corner_points.append(Point(xwidth, ywidth, 0.0))
56 corner_points.append(Point(0.0, ywidth, 0.0))
57 corner_points.append(corner_points[0]) #repeated point for line looping
58
59 # Edges of Free Surface
60 hor_lines=[]
61 for i in range(0,4): # loop through four sides
62 hor_lines.append(Line(corner_points[i],corner_points[i+1]))
63
64 # Create Free Surface
65 free_surf = PlaneSurface(CurveLoop(*tuple(hor_lines[0:4])))
66
67 # Return Surface and primative arrays.
68 return free_surf,hor_lines,corner_points
69
70 def buildLayer(xwidth,ywidth,depth,lay_surf,hor_lines,corner_points):
71 '''
72 Builds a boxlike volume and returns primatives for layered model
73 construction.
74 Parameters:
75 xwidth :: width in x direction in meters.
76 ywidth :: width in y direction in meters.
77 depth :: depth to bottom of layer in meters.
78 lay_surf :: surface at top of layer (see buildFreeSurf)
79 hor_lines :: lines of lay_surf
80 corner_points :: points of hor_lines
81 '''
82 # Layer Corners
83 corner_points.append(Point(0.0, 0.0, depth))
84 corner_points.append(Point(xwidth, 0.0, depth))
85 corner_points.append(Point(xwidth, ywidth, depth))
86 corner_points.append(Point(0.0, ywidth, depth))
87 corner_points.append(corner_points[5]) #repeated point for line looping
88
89 # Build the bottom surface edges.
90 for i in range(0,4): # loop through four edges
91 hor_lines.append(Line(corner_points[5+i],corner_points[6+i]))
92
93 # Join corners vertically.
94 ver_lines=[]
95 for i in range(0,4): # loop through four corners
96 ver_lines.append(Line(corner_points[i],corner_points[i+5]))
97 ver_lines.append(ver_lines[0]) #repeated edge for surface looping
98
99 # Build surface array.
100 lay_surf=[-lay_surf] #Negative of top surface
101 # Bottom Surface
102 lay_surf.append(PlaneSurface(CurveLoop(*tuple(hor_lines[4:8]))))
103 for i in range(0,4): # loop through four sides
104 lay_surf.append(PlaneSurface(CurveLoop(\
105 -ver_lines[i],-hor_lines[i+4],\
106 ver_lines[i+1],hor_lines[i] )))
107
108 # Build Layer Volume
109 lay_vol=Volume(-SurfaceLoop(*tuple(lay_surf)))
110
111 # Return new volume, and primatives for next volume layer.
112 return lay_vol,-lay_surf[1],hor_lines[4:8],corner_points[5:10]
113
114
115 def layer_cake(xwidth,ywidth,depths,ele_size):
116 '''
117 Builds a horizontally layered box like model. All layers are
118 tagged as 'interface_i' where i is the python style integer denoting
119 that layer. For example, the free surface is tagged 'interface_0'.
120 Volumes are similarly tagged as 'volume_i'.
121
122 Parameters:
123 xwidth :: width in x direction in meters.
124 ywidth :: width in y direction in meters.
125 depth :: depth to bottom of layer in meters.
126 ele_size :: the element meshing size.
127 fname :: the output file name.
128 KW Arg:
129 save_path :: path to save outputs.
130
131 One may save the domain using:
132 # Output settings.
133 domain.setScriptFileName(os.path.join(save_path,fname+".geo"))
134 domain.setMeshFileName(os.path.join(save_path,fname+".msh"))
135 findomain=fin.MakeDomain(domain) # make the finley domain:
136
137 # Create a file that can be read back in to python with
138 # ReadMesh
139 findomain.write(os.path.join(save_path,fname+".fly"))
140 '''
141
142 #get number of layers
143 ndepths=len(depths)
144
145 # Specify the domain.
146 domain=Design(dim=3,element_size=ele_size)
147
148 # Build the First Surface and add it to the domain
149 fsuf,fsurl,fsurp=buildFreeSurface(xwidth,ywidth)
150 domain.addItems(PropertySet('intface_%d'%(0),fsuf))
151
152 # Build each layer sequentially
153 # Set up temporary variables.
154 tsuf=fsuf; tsurl=fsurl; tsurp=fsurp
155 # Loop through all depths.
156 for i in range(0,ndepths):
157 # Build new layer.
158 tvol,tsuf,tsurl,tsurp=buildLayer(xwidth,ywidth,depths[i],\
159 tsuf,tsurl,tsurp)
160 # Add the new interface to the domain.
161 domain.addItems(PropertySet('intface_%d'%(i+1),tsuf))
162 # Add the new volume/layer to the domain.
163 domain.addItems(PropertySet('volume_%d'%i,tvol))
164
165 return domain
166

  ViewVC Help
Powered by ViewVC 1.1.26