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

Annotation of /trunk/doc/examples/cookbook/example10c_0.py

Parent Directory Parent Directory | Revision Log Revision Log


Revision 3195 - (hide annotations)
Wed Sep 22 00:28:04 2010 UTC (8 years, 4 months ago) by ahallam
Original Path: trunk/doc/examples/cookbook/faultgen.py
File MIME type: text/x-python
File size: 7695 byte(s)
Shortened runtime of cookbook examples to aid testing.
1 ahallam 3135
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 ahallam 3195 fault_dip_front=30*DEG
40     fault_dip_back=30*DEG
41 ahallam 3135
42 ahallam 3195 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 ahallam 3135 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 ahallam 3195 des=Design(dim=3, order=1, element_size = 40*m, keep_files=True)
179 ahallam 3135 des.addItems(*tuple(PS.values()))
180     des.addItems(PropertySet("fault",Volume(SurfaceLoop( *tuple(FF)))))
181 ahallam 3195 des.setMeshFileName("fault.msh")
182 ahallam 3135 dom=MakeDomain(des)
183 ahallam 3195 dom.write("fault.fly")
184 ahallam 3135

  ViewVC Help
Powered by ViewVC 1.1.26