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

  ViewVC Help
Powered by ViewVC 1.1.26