/[escript]/trunk/doc/examples/cookbook/layer_cake.py
ViewVC logotype

Annotation of /trunk/doc/examples/cookbook/layer_cake.py

Parent Directory Parent Directory | Revision Log Revision Log


Revision 3058 - (hide annotations)
Tue Jul 6 02:13:14 2010 UTC (10 years, 6 months ago) by ahallam
File MIME type: text/x-python
File size: 4711 byte(s)


1 ahallam 3058 # -*- 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     # Build a free surface for start of layer cake model.
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])
59     hor_lines=[]
60     for i in range(0,4): #top of layer loop
61     hor_lines.append(Line(corner_points[i],corner_points[i+1]))
62     #free surface
63     free_surf = PlaneSurface(CurveLoop(*tuple(hor_lines[0:4])))
64     return free_surf,hor_lines,corner_points
65    
66     def buildLayer(xwidth,ywidth,depth,lay_surf,hor_lines,corner_points):
67     # Layer Corners
68     corner_points.append(Point(0.0, 0.0, depth))
69     corner_points.append(Point(xwidth, 0.0, depth))
70     corner_points.append(Point(xwidth, ywidth, depth))
71     corner_points.append(Point(0.0, ywidth, depth))
72     corner_points.append(corner_points[5])
73     print corner_points
74    
75     # Join corners in anti-clockwise manner.
76     # for i in range(0,4): #top of layer loop
77     # hor_lines.append(Line(corner_points[i],corner_points[i+1])
78     for i in range(0,4): #bottom of layer loop
79     hor_lines.append(Line(corner_points[5+i],corner_points[6+i]))
80     # Join corners vertically.
81     ver_lines=[]
82     for i in range(0,5):
83     ver_lines.append(Line(corner_points[i],corner_points[i+5]))
84     print 'corners',i,i+5
85    
86     #bottom of layer
87     lay_surf=[PlaneSurface(CurveLoop(-hor_lines[3],-hor_lines[2],\
88     -hor_lines[1],-hor_lines[1]))]
89     lay_surf.append(PlaneSurface(CurveLoop(*tuple(hor_lines[4:8]))))
90     # lay_subot=PlaneSurface(-CurveLoop(*tuple(hor_lines[4:8])))
91     #sides of layer
92     for i in range(0,4):
93     lay_surf.append(PlaneSurface(CurveLoop(\
94     -ver_lines[i],\
95     -hor_lines[i+4],\
96     ver_lines[i+1],\
97     hor_lines[i] )))
98     lay_vol=Volume(-SurfaceLoop(*tuple(lay_surf)))
99     #lay_vol=lay_surf
100     return lay_vol,lay_surf[1],hor_lines[4:8],corner_points[5:10]
101    
102     def LayerCake(xwidth,ywidth,depths,ele_size,fname,save_path=""):
103     #get number of layers
104     ndepths=len(depths)
105     fsuf,fsurl,fsurp=buildFreeSurface(xwidth,ywidth)
106     domain=Design(dim=3,element_size=ele_size)
107     #domain.addItems(PropertySet('free_surf',fsuf))
108     #build each layer sequentially
109     tsuf=fsuf; tsurl=fsurl; tsurp=fsurp
110     for i in range(0,ndepths):
111     tvol,tsuf,tsurl,tsurp=buildLayer(xwidth,ywidth,depths[i],\
112     tsuf,tsurl,tsurp)
113     #domain.addItems(*tuple(tvol))
114     #domain.addItems(PropertySet('intface_%d'%(i+1),tsuf))
115     domain.addItems(PropertySet('volume_%d'%i,tvol))
116     # Add the subdomains and flux boundaries.
117    
118     domain.setScriptFileName(os.path.join(save_path,fname+".geo"))
119     domain.setMeshFileName(os.path.join(save_path,fname+".msh"))
120     findomain=MakeDomain(domain) # make the finley domain:
121     # Create a file that can be read back in to python with
122     #mesh=ReadMesh(fileName)
123     #findomain.write(os.path.join(save_path,fname+".fly"))
124    
125     LayerCake(100.0,100.0,[100.],10.,'testlc')
126    

  ViewVC Help
Powered by ViewVC 1.1.26