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

  ViewVC Help
Powered by ViewVC 1.1.26