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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 3060 - (show annotations)
Tue Jul 13 00:55:09 2010 UTC (9 years, 4 months ago) by ahallam
File MIME type: text/x-python
File size: 6202 byte(s)
Adding layer cake builder to 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 from esys.finley import MakeDomain #Converter for escript
42 from esys.escript.unitsSI import *
43 import os
44 from math import *
45 import pylab as pl
46 import numpy as np
47
48
49
50 def buildFreeSurface(xwidth,ywidth):
51 '''
52 Build a free surface to start the layer cake model.
53 This surface is planar.
54 Parameters:
55 xwidth :: width in x direction in meters.
56 ywidth :: width in y direction in meters.
57 '''
58
59 # Layer Corners
60 corner_points=[]
61 corner_points.append(Point(0.0, 0.0, 0.0))
62 corner_points.append(Point(xwidth, 0.0, 0.0))
63 corner_points.append(Point(xwidth, ywidth, 0.0))
64 corner_points.append(Point(0.0, ywidth, 0.0))
65 corner_points.append(corner_points[0]) #repeated point for line looping
66
67 # Edges of Free Surface
68 hor_lines=[]
69 for i in range(0,4): # loop through four sides
70 hor_lines.append(Line(corner_points[i],corner_points[i+1]))
71
72 # Create Free Surface
73 free_surf = PlaneSurface(CurveLoop(*tuple(hor_lines[0:4])))
74
75 # Return Surface and primative arrays.
76 return free_surf,hor_lines,corner_points
77
78 def buildLayer(xwidth,ywidth,depth,lay_surf,hor_lines,corner_points):
79 '''
80 Builds a boxlike volume and returns primatives for layered model
81 construction.
82 Parameters:
83 xwidth :: width in x direction in meters.
84 ywidth :: width in y direction in meters.
85 depth :: depth to bottom of layer in meters.
86 lay_surf :: surface at top of layer (see buildFreeSurf)
87 hor_lines :: lines of lay_surf
88 corner_points :: points of hor_lines
89 '''
90 # Layer Corners
91 corner_points.append(Point(0.0, 0.0, depth))
92 corner_points.append(Point(xwidth, 0.0, depth))
93 corner_points.append(Point(xwidth, ywidth, depth))
94 corner_points.append(Point(0.0, ywidth, depth))
95 corner_points.append(corner_points[5]) #repeated point for line looping
96
97 # Build the bottom surface edges.
98 for i in range(0,4): # loop through four edges
99 hor_lines.append(Line(corner_points[5+i],corner_points[6+i]))
100
101 # Join corners vertically.
102 ver_lines=[]
103 for i in range(0,4): # loop through four corners
104 ver_lines.append(Line(corner_points[i],corner_points[i+5]))
105 ver_lines.append(ver_lines[0]) #repeated edge for surface looping
106
107 # Build surface array.
108 lay_surf=[-lay_surf] #Negative of top surface
109 # Bottom Surface
110 lay_surf.append(PlaneSurface(CurveLoop(*tuple(hor_lines[4:8]))))
111 for i in range(0,4): # loop through four sides
112 lay_surf.append(PlaneSurface(CurveLoop(\
113 -ver_lines[i],-hor_lines[i+4],\
114 ver_lines[i+1],hor_lines[i] )))
115
116 # Build Layer Volume
117 lay_vol=Volume(-SurfaceLoop(*tuple(lay_surf)))
118
119 # Return new volume, and primatives for next volume layer.
120 return lay_vol,-lay_surf[1],hor_lines[4:8],corner_points[5:10]
121
122
123 def LayerCake(xwidth,ywidth,depths,ele_size,fname,save_path=""):
124 '''
125 Builds a horizontally layered box like model. All layers are
126 tagged as 'interface_i' where i is the python style integer denoting
127 that layer. For example, the free surface is tagged 'interface_0'.
128 Volumes are similarly tagged as 'volume_i'.
129
130 Parameters:
131 xwidth :: width in x direction in meters.
132 ywidth :: width in y direction in meters.
133 depth :: depth to bottom of layer in meters.
134 ele_size :: the element meshing size.
135 fname :: the output file name.
136 KW Arg:
137 save_path :: path to save outputs.
138 '''
139
140 #get number of layers
141 ndepths=len(depths)
142
143 # Specify the domain.
144 domain=Design(dim=3,element_size=ele_size)
145
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 # Output settings.
164 domain.setScriptFileName(os.path.join(save_path,fname+".geo"))
165 domain.setMeshFileName(os.path.join(save_path,fname+".msh"))
166 findomain=MakeDomain(domain) # make the finley domain:
167
168 # Create a file that can be read back in to python with
169 # ReadMesh
170 findomain.write(os.path.join(save_path,fname+".fly"))
171
172 #LayerCake(100.0,100.0,[10.,40.,80.,100.,150.],10.,'testlc')
173

  ViewVC Help
Powered by ViewVC 1.1.26