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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 3069 - (hide annotations)
Wed Jul 21 03:24:48 2010 UTC (11 years, 3 months ago) by ahallam
File MIME type: text/x-python
File size: 5930 byte(s)
Updates to layer cake and cookbook as well as examples.
1 ahallam 3060 # -*- 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 ahallam 3069
47     Parameters:
48     xwidth :: width in x direction in meters.\
49     ywidth :: width in y direction in meters.\
50 ahallam 3060 '''
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 ahallam 3069
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 ahallam 3060 '''
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 ahallam 3069 def layer_cake(domain,xwidth,ywidth,depths):
118 ahallam 3060 '''
119     Builds a horizontally layered box like model. All layers are
120     tagged as 'interface_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 ahallam 3069 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 ahallam 3061
130     One may save the domain using:
131 ahallam 3069 # 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 ahallam 3061
136 ahallam 3069 Create a file that can be read back in to python with ReadMesh.\
137     findomain.write(os.path.join(save_path,fname+".fly"))\
138 ahallam 3060 '''
139    
140     #get number of layers
141     ndepths=len(depths)
142    
143 ahallam 3069 if domain.getDim() <> 3:
144     raise TypeError("domain must be of dimension order 3.")
145 ahallam 3060
146     # Build the First Surface and add it to the domain
147     fsuf,fsurl,fsurp=buildFreeSurface(xwidth,ywidth)
148     domain.addItems(PropertySet('intface_%d'%(0),fsuf))
149    
150     # Build each layer sequentially
151     # Set up temporary variables.
152     tsuf=fsuf; tsurl=fsurl; tsurp=fsurp
153     # Loop through all depths.
154     for i in range(0,ndepths):
155     # Build new layer.
156     tvol,tsuf,tsurl,tsurp=buildLayer(xwidth,ywidth,depths[i],\
157     tsuf,tsurl,tsurp)
158     # Add the new interface to the domain.
159     domain.addItems(PropertySet('intface_%d'%(i+1),tsuf))
160     # Add the new volume/layer to the domain.
161     domain.addItems(PropertySet('volume_%d'%i,tvol))
162    
163 ahallam 3061 return domain
164 ahallam 3060

  ViewVC Help
Powered by ViewVC 1.1.26