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

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

Parent Directory Parent Directory | Revision Log Revision Log


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


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 # 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