/[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 3774 - (show annotations)
Wed Jan 18 06:29:34 2012 UTC (7 years, 9 months ago) by jfenwick
File MIME type: text/x-python
File size: 5987 byte(s)
dudley, pasowrap, pycad

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
47 Parameters:
48 xwidth :: width in x direction in meters.\
49 ywidth :: width in y direction in meters.\
50 '''
51
52 # Layer Corners
53 corner_points=[]
54 corner_points.append(Point(0.0, 0.0, 0.0))
55 corner_points.append(Point(xwidth, 0.0, 0.0))
56 corner_points.append(Point(xwidth, ywidth, 0.0))
57 corner_points.append(Point(0.0, ywidth, 0.0))
58 corner_points.append(corner_points[0]) #repeated point for line looping
59
60 # Edges of Free Surface
61 hor_lines=[]
62 for i in range(0,4): # loop through four sides
63 hor_lines.append(Line(corner_points[i],corner_points[i+1]))
64
65 # Create Free Surface
66 free_surf = PlaneSurface(CurveLoop(*tuple(hor_lines[0:4])))
67
68 # Return Surface and primative arrays.
69 return free_surf,hor_lines,corner_points
70
71 def buildLayer(xwidth,ywidth,depth,lay_surf,hor_lines,corner_points):
72 '''
73 Builds a boxlike volume and returns primatives for layered model
74 construction.
75
76 Parameters:
77 xwidth :: width in x direction in meters.\
78 ywidth :: width in y direction in meters.\
79 depth :: depth to bottom of layer in meters.\
80 lay_surf :: surface at top of layer (see buildFreeSurf)\
81 hor_lines :: lines of lay_surf\
82 corner_points :: points of hor_lines\
83 '''
84 # Layer Corners
85 corner_points.append(Point(0.0, 0.0, depth))
86 corner_points.append(Point(xwidth, 0.0, depth))
87 corner_points.append(Point(xwidth, ywidth, depth))
88 corner_points.append(Point(0.0, ywidth, depth))
89 corner_points.append(corner_points[5]) #repeated point for line looping
90
91 # Build the bottom surface edges.
92 for i in range(0,4): # loop through four edges
93 hor_lines.append(Line(corner_points[5+i],corner_points[6+i]))
94
95 # Join corners vertically.
96 ver_lines=[]
97 for i in range(0,4): # loop through four corners
98 ver_lines.append(Line(corner_points[i],corner_points[i+5]))
99 ver_lines.append(ver_lines[0]) #repeated edge for surface looping
100
101 # Build surface array.
102 lay_surf=[-lay_surf] #Negative of top surface
103 # Bottom Surface
104 lay_surf.append(PlaneSurface(CurveLoop(*tuple(hor_lines[4:8]))))
105 for i in range(0,4): # loop through four sides
106 lay_surf.append(PlaneSurface(CurveLoop(\
107 -ver_lines[i],-hor_lines[i+4],\
108 ver_lines[i+1],hor_lines[i] )))
109
110 # Build Layer Volume
111 lay_vol=Volume(-SurfaceLoop(*tuple(lay_surf)))
112
113 # Return new volume, and primatives for next volume layer.
114 return lay_vol,-lay_surf[1],hor_lines[4:8],corner_points[5:10]
115
116
117 def layer_cake(domain,xwidth,ywidth,depths):
118 '''
119 Builds a horizontally layered box like model. All layers are
120 tagged as 'intface_i' where i is the python style integer denoting
121 that layer. For example, the free surface is tagged 'interface_0'.
122 Volumes are similarly tagged as 'volume_i'.
123
124 Parameters:
125 domain :: output of Pycad.Design - It needs to be dim 3.\
126 xwidth :: width in x direction in meters.\
127 ywidth :: width in y direction in meters.\
128 depth :: depth to bottom of layer in meters.\
129
130 One may save the domain using:
131 # Output settings.\
132 domain.setScriptFileName(os.path.join(save_path,fname+".geo"))\
133 domain.setMeshFileName(os.path.join(save_path,fname+".msh"))\
134 findomain=fin.MakeDomain(domain) # make the finley domain:\
135
136 Create a file that can be read back in to python with ReadMesh.\
137 findomain.write(os.path.join(save_path,fname+".fly"))\
138 '''
139
140 #get number of layers
141 if not hasattr(depths,"__len__"): depths = [ depths, ]
142 ndepths=len(depths)
143
144 if domain.getDim() != 3:
145 raise TypeError("domain must be of dimension order 3.")
146
147 # Build the First Surface and add it to the domain
148 fsuf,fsurl,fsurp=buildFreeSurface(xwidth,ywidth)
149 domain.addItems(PropertySet('intface_%d'%(0),fsuf))
150
151 # Build each layer sequentially
152 # Set up temporary variables.
153 tsuf=fsuf; tsurl=fsurl; tsurp=fsurp
154 # Loop through all depths.
155 for i in range(0,ndepths):
156 # Build new layer.
157 tvol,tsuf,tsurl,tsurp=buildLayer(xwidth,ywidth,depths[i],\
158 tsuf,tsurl,tsurp)
159 # Add the new interface to the domain.
160 domain.addItems(PropertySet('intface_%d'%(i+1),tsuf))
161 # Add the new volume/layer to the domain.
162 domain.addItems(PropertySet('volume_%d'%i,tvol))
163
164 return domain
165

  ViewVC Help
Powered by ViewVC 1.1.26