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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 5707 - (show annotations)
Mon Jun 29 03:59:06 2015 UTC (3 years, 8 months ago) by sshaw
File MIME type: text/x-python
File size: 6403 byte(s)
adding copyright headers to files without copyright info, moved header to top of file in some cases where it wasn't
1 ##############################################################################
2 #
3 # Copyright (c) 2009-2015 by The University of Queensland
4 # http://www.uq.edu.au
5 #
6 # Primary Business: Queensland, Australia
7 # Licensed under the Open Software License version 3.0
8 # http://www.opensource.org/licenses/osl-3.0.php
9 #
10 # Development until 2012 by Earth Systems Science Computational Center (ESSCC)
11 # Development 2012-2013 by School of Earth Sciences
12 # Development from 2014 by Centre for Geoscience Computing (GeoComp)
13 #
14 ##############################################################################
15 from __future__ import division, print_function
16
17 __copyright__="""Copyright (c) 2009-2015 by The 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 Author: Antony Hallam antony.hallam@uqconnect.edu.au
26 """
27
28
29 ############################################################FILE HEADER
30 # example09n.py
31 # Create a simple 3D model for use in example09. This is the low res
32 # mesh for illustration purposes only.
33 #
34 #######################################################EXTERNAL MODULES
35 from esys.pycad import * #domain constructor
36 from esys.pycad.gmsh import Design #Finite Element meshing package
37 from esys.escript import mkDir, getMPISizeWorld
38 import os
39 import math
40 try:
41 # This imports the rectangle domain function
42 from esys.finley import MakeDomain
43 HAVE_FINLEY = True
44 except ImportError:
45 print("Finley module not available")
46 HAVE_FINLEY = False
47 ########################################################MPI WORLD CHECK
48 if getMPISizeWorld() > 1:
49 import sys
50 print("This example will not run in an MPI world.")
51 sys.exit(0)
52
53 if HAVE_FINLEY:
54 # make sure path exists
55 save_path= os.path.join("data","example09c")
56 mkDir(save_path)
57
58 ################################################ESTABLISHING PARAMETERS
59 #Model Parameters
60 origin=[0,0] #orign of model
61 xwidth=300.0 #width of model
62 depth=-100.0 #depth of model
63 nintf=3 #number of the interfaces
64 lintf_depths=[-20,-40,-60] #depth of interfaces
65 rintf_depths=[-30,-50,-70] #vertical displacement across fault
66 fault_dip=40.0 #dip of fault plane
67 fault_atsurface=50.0 #location of fault at surface
68 fault_width=10 #apparent width of fault plane
69
70 element_size=1.
71
72 ####################################################DOMAIN CONSTRUCTION
73 x0=0.0+origin[0]
74 y0=0.0+origin[1]
75 #z=0.0+origin[2]
76 fault_atsurface=fault_atsurface+origin[0]
77 xwidth=xwidth+origin[0]
78 depth=depth+origin[1]
79 # Construction points, 4 vectors that descend from the surface with nintf+2 points.
80 left_edge=[Point(x0,y0)];
81 leftf_edge=[Point(fault_atsurface,y0)];
82 rightf_edge=[Point(fault_atsurface+fault_width,y0)];
83 right_edge=[Point(xwidth,y0)];
84
85 for i in range(0,nintf):
86 left_shift=(lintf_depths[i]-y0)/math.tan(fault_dip)
87 right_shift=(rintf_depths[i]-y0)/math.tan(fault_dip)
88 left_edge.append(Point(x0,lintf_depths[i]+origin[1]))
89 leftf_edge.append(Point(fault_atsurface+left_shift,lintf_depths[i]+origin[1]))
90 rightf_edge.append(Point(fault_atsurface+right_shift+fault_width,rintf_depths[i]+origin[1]))
91 right_edge.append(Point(xwidth,rintf_depths[i]+origin[1]))
92
93 left_edge.append(Point(x0,depth))
94 leftf_edge.append(Point(fault_atsurface+(depth-y0)/math.tan(fault_dip),depth))
95 rightf_edge.append(Point(fault_atsurface+fault_width+(depth-y0)/math.tan(fault_dip),depth))
96 right_edge.append(Point(xwidth,depth))
97
98 #Build lines
99 lright=[]; nlright=[];
100 lhright=[]; nlhright=[];
101 lfright=[]; nlfright=[];
102 lfhor=[]; nlfhor=[];
103 lfleft=[]; nlfleft=[];
104 lhleft=[]; nlhleft=[];
105 lleft=[]; nlleft=[];
106
107 #Build vertical lines
108 for i in range(0,nintf+1):
109 lleft.append(Line(left_edge[i],left_edge[i+1]))
110 lfleft.append(Line(leftf_edge[i],leftf_edge[i+1]))
111 lfright.append(Line(rightf_edge[i],rightf_edge[i+1]))
112 lright.append(Line(right_edge[i],right_edge[i+1]))
113
114 #Build horizontal lines
115 for i in range(0,nintf+2):
116 lhleft.append(Line(left_edge[i],leftf_edge[i]))
117 lhright.append(Line(rightf_edge[i],right_edge[i]))
118 lfhor.append(Line(leftf_edge[0],rightf_edge[0]))
119 lfhor.append(Line(leftf_edge[nintf+1],rightf_edge[nintf+1]))
120
121 #Build negative lines
122 for i in range(nintf,-1,-1):
123 nlleft.append(-lleft[i])
124 nlfleft.append(-lfleft[i])
125 nlfright.append(-lfright[i])
126 nlright.append(-lright[i])
127 for i in range(nintf+1,-1,-1):
128 nlhleft.append(-lhleft[i])
129 nlhright.append(-lhright[i])
130
131 #Build curveloops
132 lcurves=[]
133 fcurves=[]
134 rcurves=[]
135
136 #Fault
137 for i in range(0,nintf+1):
138 fcurves.append(lfleft[i])
139 fcurves.append(lfhor[1])
140 for i in range(0,nintf+1):
141 fcurves.append(nlfright[i])
142 fcurves.append(-lfhor[0])
143 fcurves=CurveLoop(*tuple(fcurves))
144
145 #Left and Right Blocks
146 for i in range(0,nintf+1):
147 lcurves.append(CurveLoop(lleft[i],lhleft[i+1],nlfleft[nintf-i],nlhleft[nintf+1-i]))
148 rcurves.append(CurveLoop(lfright[i],lhright[i+1],nlright[nintf-i],nlhright[nintf+1-i]))
149
150 #Build Surfaces
151 fsurf=PlaneSurface(fcurves)
152 lsurf=[]
153 rsurf=[]
154 for i in range(0,nintf+1):
155 lsurf.append(PlaneSurface(lcurves[i]))
156 rsurf.append(PlaneSurface(rcurves[i]))
157
158 d=Design(dim=2, element_size=element_size, order=2)
159
160 d.addItems(PropertySet('fault',fsurf))
161 for i in range(0,nintf+1):
162 d.addItems(PropertySet('lblock%d'%i,lsurf[i]))
163 d.addItems(PropertySet('rblock%d'%i,rsurf[i]))
164
165 d.addItems(PropertySet('top',lhright[0],lfhor[0],lhleft[0],lhright[4],lhleft[4],lfhor[1]))
166
167 d.setScriptFileName(os.path.join(save_path,"example09n.geo"))
168 d.setMeshFileName(os.path.join(save_path,"example09n.msh"))
169 #
170 # make the domain:
171 #
172 domain=MakeDomain(d)
173 # mesh=ReadMesh(fileName) this is how to read the fly file into escript
174 domain.write(os.path.join(save_path,"example09n.fly"))
175
176
177
178
179

  ViewVC Help
Powered by ViewVC 1.1.26