/[escript]/branches/mvexpr/pycad/py_src/extras.py
ViewVC logotype

Contents of /branches/mvexpr/pycad/py_src/extras.py

Parent Directory Parent Directory | Revision Log Revision Log


Revision 4401 - (show annotations)
Fri May 10 05:14:06 2013 UTC (6 years, 5 months ago) by jfenwick
File MIME type: text/x-python
File size: 6049 byte(s)
Trying a directory swap

1 # -*- coding: utf-8 -*-
2
3 ##############################################################################
4 #
5 # Copyright (c) 2003-2013 by University of Queensland
6 # http://www.uq.edu.au
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 # Development until 2012 by Earth Systems Science Computational Center (ESSCC)
13 # Development since 2012 by School of Earth Sciences
14 #
15 ##############################################################################
16
17 __copyright__="""Copyright (c) 2003-2013 by 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 Layer Cake
26
27 This script uses the pycad module to build a rectangular layered
28 system of prisms for modelling.
29
30 :var __author__: name of author
31 :var __copyright__: copyrights
32 :var __license__: licence agreement
33 :var __url__: url entry point on documentation
34 :var __version__: version
35 :var __date__: date of the version
36 """
37
38 __author__="Antony Hallam antony.hallam@uqconnect.edu.au"
39
40 from esys.pycad import * #domain constructor
41 from esys.pycad.gmsh import Design #Finite Element meshing package
42
43 def buildFreeSurface(xwidth,ywidth):
44 '''
45 Build a free surface to start the layer cake model.
46 This surface is planar.
47
48 Parameters:
49 xwidth :: width in x direction in meters.\
50 ywidth :: width in y direction in meters.\
51 '''
52
53 # Layer Corners
54 corner_points=[]
55 corner_points.append(Point(0.0, 0.0, 0.0))
56 corner_points.append(Point(xwidth, 0.0, 0.0))
57 corner_points.append(Point(xwidth, ywidth, 0.0))
58 corner_points.append(Point(0.0, ywidth, 0.0))
59 corner_points.append(corner_points[0]) #repeated point for line looping
60
61 # Edges of Free Surface
62 hor_lines=[]
63 for i in range(0,4): # loop through four sides
64 hor_lines.append(Line(corner_points[i],corner_points[i+1]))
65
66 # Create Free Surface
67 free_surf = PlaneSurface(CurveLoop(*tuple(hor_lines[0:4])))
68
69 # Return Surface and primative arrays.
70 return free_surf,hor_lines,corner_points
71
72 def buildLayer(xwidth,ywidth,depth,lay_surf,hor_lines,corner_points):
73 '''
74 Builds a boxlike volume and returns primatives for layered model
75 construction.
76
77 Parameters:
78 xwidth :: width in x direction in meters.\
79 ywidth :: width in y direction in meters.\
80 depth :: depth to bottom of layer in meters.\
81 lay_surf :: surface at top of layer (see buildFreeSurf)\
82 hor_lines :: lines of lay_surf\
83 corner_points :: points of hor_lines\
84 '''
85 # Layer Corners
86 corner_points.append(Point(0.0, 0.0, depth))
87 corner_points.append(Point(xwidth, 0.0, depth))
88 corner_points.append(Point(xwidth, ywidth, depth))
89 corner_points.append(Point(0.0, ywidth, depth))
90 corner_points.append(corner_points[5]) #repeated point for line looping
91
92 # Build the bottom surface edges.
93 for i in range(0,4): # loop through four edges
94 hor_lines.append(Line(corner_points[5+i],corner_points[6+i]))
95
96 # Join corners vertically.
97 ver_lines=[]
98 for i in range(0,4): # loop through four corners
99 ver_lines.append(Line(corner_points[i],corner_points[i+5]))
100 ver_lines.append(ver_lines[0]) #repeated edge for surface looping
101
102 # Build surface array.
103 lay_surf=[-lay_surf] #Negative of top surface
104 # Bottom Surface
105 lay_surf.append(PlaneSurface(CurveLoop(*tuple(hor_lines[4:8]))))
106 for i in range(0,4): # loop through four sides
107 lay_surf.append(PlaneSurface(CurveLoop(\
108 -ver_lines[i],-hor_lines[i+4],\
109 ver_lines[i+1],hor_lines[i] )))
110
111 # Build Layer Volume
112 lay_vol=Volume(-SurfaceLoop(*tuple(lay_surf)))
113
114 # Return new volume, and primatives for next volume layer.
115 return lay_vol,-lay_surf[1],hor_lines[4:8],corner_points[5:10]
116
117
118 def layer_cake(domain,xwidth,ywidth,depths):
119 '''
120 Builds a horizontally layered box like model. All layers are
121 tagged as 'intface_i' where i is the python style integer denoting
122 that layer. For example, the free surface is tagged 'interface_0'.
123 Volumes are similarly tagged as 'volume_i'.
124
125 Parameters:
126 domain :: output of Pycad.Design - It needs to be dim 3.\
127 xwidth :: width in x direction in meters.\
128 ywidth :: width in y direction in meters.\
129 depth :: depth to bottom of layer in meters.\
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 ReadMesh.\
138 findomain.write(os.path.join(save_path,fname+".fly"))\
139 '''
140
141 #get number of layers
142 if not hasattr(depths,"__len__"): depths = [ depths, ]
143 ndepths=len(depths)
144
145 if domain.getDim() != 3:
146 raise TypeError("domain must be of dimension order 3.")
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