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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 3195 - (show annotations)
Wed Sep 22 00:28:04 2010 UTC (9 years, 2 months ago) by ahallam
File MIME type: text/x-python
File size: 7695 byte(s)
Shortened runtime of cookbook examples to aid testing.
1
2 ########################################################
3 #
4 # Copyright (c) 2003-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) 2003-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 from esys.escript import *
23 from esys.pycad import *
24 from esys.pycad.gmsh import Design
25 from esys.finley import MakeDomain
26 from esys.escript.unitsSI import *
27 from math import tan
28
29 #
30 # input data:
31 #
32 width=5*km
33 length=5*km
34 depth=2*km
35 fault_w_front=200*m
36 fault_w_back=200*m
37 fault_mid_front=2*km
38 fault_mid_back=4*km
39 fault_dip_front=30*DEG
40 fault_dip_back=30*DEG
41
42 layers_left_at_front=[ [ 'limestone' , 1100*m ] , [ 'xx' , 150*m ] ,['shale',300*m], [ 'limestone' ] ]
43 layers_right_at_front=[ [ 'limestone' , 500*m ], [ 'xx' , 150*m ],['shale',300*m], ['limestone' ] ]
44 slop=2*DEG
45 #
46 # ====================================================================================
47 #
48 def setOffsetsFromThickness(layers):
49 l=0
50 offsets=[]
51 for n in layers:
52 if len(n)<2:
53 d=depth-l
54 if depth-l <=0:
55 raise "Layer structure is too thick. Increase depth > %e"%l
56 else:
57 d=n[1]
58 l+=d
59 offsets.append([ n[0], -(l-d), d, -l ])
60 return offsets
61 def setVerticalPositionsFromDip(layers, dip):
62 offsets=[]
63 s=1/tan(90*DEG-dip)
64 for n in layers:
65 name=n[0]
66 top=n[1]
67 d=n[2]
68 bot=n[3]
69 offsets.append([ name, top, s*top, d, bot, s*bot ] )
70 return offsets
71
72 def setBackLayers(layers, slop, dip):
73 if layers[0][3]-tan(slop)*length <=0:
74 raise ValueError,"negative thickness in %s at back"%(layers[0][0],)
75 out=[[ layers[0][0], layers[0][3]-tan(slop)*length ]]
76 for i in range(1,len(layers)-1): out.append([layers[i][0], layers[i][3] ])
77 out.append([ layers[-1][0] ])
78 out=setOffsetsFromThickness(out)
79 return setVerticalPositionsFromDip(out, dip)
80
81 def getCutLine(p0,layers,offset=True):
82 out=[ ]
83 p=p0
84 for n in layers:
85 if offset:
86 p1,p = p, Point(p0.getCoordinates()[0]+n[5],p0.getCoordinates()[1],p0.getCoordinates()[2]+n[4])
87 else:
88 p1,p = p, Point(p0.getCoordinates()[0],p0.getCoordinates()[1],p0.getCoordinates()[2]+n[4])
89 out.append(Line(p1,p))
90 return out
91
92 def getBackLine(front, back):
93 out = [ Line(front[0].getStartPoint(), back[0].getStartPoint()) ]
94 for i in range(len(front)):
95 out.append(Line(front[i].getEndPoint(), back[i].getEndPoint()))
96 return out
97
98 def getCrossLine(left, right):
99 return getBackLine(left,right)
100
101 def addVolume(front_left, back_left, front_right, back_right, PS, FF, map, filter_left=False):
102
103 front_to_back_left = getBackLine(front_left , back_left)
104 front_to_back_right = getBackLine(front_right , back_right)
105 front_left_to_right= getCrossLine(front_left, front_right)
106 back_left_to_right= getCrossLine(back_left, back_right)
107
108 #if filter_left:
109 # out1=front_to_back_left[0]
110 # out2=front_to_back_left[-1]
111 #else:
112 out1=front_to_back_right[0]
113 out2=front_to_back_right[-1]
114
115 topface=PlaneSurface(CurveLoop(front_left_to_right[0], front_to_back_right[0], -back_left_to_right[0], -front_to_back_left[0]))
116 for i in range(len(map)):
117 name=map[i][0]
118 face2=PlaneSurface(CurveLoop(front_left_to_right[i+1], front_to_back_right[i+1], -back_left_to_right[i+1], -front_to_back_left[i+1]))
119 face_front=PlaneSurface(CurveLoop(front_left_to_right[i+1], -front_right[i], -front_left_to_right[i], front_left[i]))
120 face_back=PlaneSurface(CurveLoop(back_left_to_right[i+1], -back_right[i], -back_left_to_right[i], back_left[i]))
121 face_left=PlaneSurface(CurveLoop(front_left[i], front_to_back_left[i+1], -back_left[i], -front_to_back_left[i]))
122 face_right=PlaneSurface(CurveLoop(front_right[i], front_to_back_right[i+1], -back_right[i], -front_to_back_right[i]))
123
124 v=Volume(SurfaceLoop(topface,-face2, face_front, -face_back, -face_left, face_right))
125 print v
126 topface=face2
127 if filter_left:
128 FF.append(face_right)
129 else:
130 FF.append(-face_right)
131 if not PS.has_key(name): PS[name]=PropertySet(name)
132 PS[name].addItem(v)
133 return out1, out2, PS, FF
134
135
136
137 layers_left_at_edge_at_front=setOffsetsFromThickness(layers_left_at_front)
138 layers_right_at_edge_at_front=setOffsetsFromThickness(layers_right_at_front)
139
140 layers_left_at_front=setVerticalPositionsFromDip(layers_left_at_edge_at_front, fault_dip_front)
141 layers_right_at_front=setVerticalPositionsFromDip(layers_right_at_edge_at_front, fault_dip_front)
142
143
144 layers_left_at_back=setBackLayers(layers_left_at_front, slop, fault_dip_back)
145 layers_right_at_back=setBackLayers(layers_right_at_front, slop, fault_dip_back)
146
147
148 left_front_edge=getCutLine(Point(0.,0.,0.), layers_left_at_front, offset=False)
149 left_front_fault=getCutLine(Point(fault_mid_front-fault_w_front/2,0.,0.), layers_left_at_front, offset=True)
150
151 right_front_fault=getCutLine(Point(fault_mid_front+fault_w_front/2,0.,0.), layers_right_at_front, offset=True)
152 right_front_edge=getCutLine(Point(width,0.,0.), layers_right_at_front, offset=False)
153
154 left_back_edge=getCutLine(Point(0.,length,0.), layers_left_at_back, offset=False)
155 left_back_fault=getCutLine(Point(fault_mid_back-fault_w_back/2,length,0.), layers_left_at_back, offset=True)
156
157 right_back_fault=getCutLine(Point(fault_mid_back+fault_w_back/2,length,0.), layers_right_at_back, offset=True)
158 right_back_edge=getCutLine(Point(width,length,0.), layers_right_at_back, offset=False)
159
160 PS={}
161 FF=[]
162 front_to_back_left_top, front_to_back_left_bot, PS, FF=addVolume(left_front_edge, left_back_edge, left_front_fault, left_back_fault, PS, FF, layers_left_at_front, filter_left=False)
163 front_to_back_right_top, front_to_back_right_bot, PS, FF=addVolume(right_front_edge, right_back_edge, right_front_fault, right_back_fault, PS, FF, layers_right_at_front, filter_left=True)
164
165 fault_line_top_front=Line(front_to_back_left_top.getStartPoint(), front_to_back_right_top.getStartPoint())
166 fault_line_bot_front=Line(front_to_back_left_bot.getStartPoint(), front_to_back_right_bot.getStartPoint())
167 fault_line_top_back=Line(front_to_back_left_top.getEndPoint(), front_to_back_right_top.getEndPoint())
168 fault_line_bot_back=Line(front_to_back_left_bot.getEndPoint(), front_to_back_right_bot.getEndPoint())
169
170 FF.append(PlaneSurface(CurveLoop(front_to_back_left_top,fault_line_top_back,-front_to_back_right_top,-fault_line_top_front)))
171 FF.append(-PlaneSurface(CurveLoop(front_to_back_left_bot,fault_line_bot_back,-front_to_back_right_bot,-fault_line_bot_front)))
172
173 FF.append(PlaneSurface(CurveLoop(*tuple([ -fault_line_top_front,fault_line_bot_front ]+left_front_fault+[ -l for l in right_front_fault ]))))
174 FF.append(-PlaneSurface(CurveLoop(*tuple([ -fault_line_top_back,fault_line_bot_back ]+left_back_fault+[ -l for l in right_back_fault ]))))
175
176
177 # war 120
178 des=Design(dim=3, order=1, element_size = 40*m, keep_files=True)
179 des.addItems(*tuple(PS.values()))
180 des.addItems(PropertySet("fault",Volume(SurfaceLoop( *tuple(FF)))))
181 des.setMeshFileName("fault.msh")
182 dom=MakeDomain(des)
183 dom.write("fault.fly")
184

  ViewVC Help
Powered by ViewVC 1.1.26