/[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 6651 - (hide annotations)
Wed Feb 7 02:12:08 2018 UTC (12 months, 2 weeks ago) by jfenwick
File MIME type: text/x-python
File size: 8331 byte(s)
Make everyone sad by touching all the files

Copyright dates update

1 jfenwick 3981 ##############################################################################
2 ahallam 3135 #
3 jfenwick 6651 # Copyright (c) 2003-2018 by The University of Queensland
4 jfenwick 3981 # http://www.uq.edu.au
5 ahallam 3135 #
6     # Primary Business: Queensland, Australia
7 jfenwick 6112 # Licensed under the Apache License, version 2.0
8     # http://www.apache.org/licenses/LICENSE-2.0
9 ahallam 3135 #
10 jfenwick 3981 # Development until 2012 by Earth Systems Science Computational Center (ESSCC)
11 jfenwick 4657 # Development 2012-2013 by School of Earth Sciences
12     # Development from 2014 by Centre for Geoscience Computing (GeoComp)
13 jfenwick 3981 #
14     ##############################################################################
15 sshaw 5707 from __future__ import division, print_function
16 ahallam 3135
17 jfenwick 6651 __copyright__="""Copyright (c) 2003-2018 by The University of Queensland
18 jfenwick 3981 http://www.uq.edu.au
19 ahallam 3135 Primary Business: Queensland, Australia"""
20 jfenwick 6112 __license__="""Licensed under the Apache License, version 2.0
21     http://www.apache.org/licenses/LICENSE-2.0"""
22 ahallam 3135 __url__="https://launchpad.net/escript-finley"
23    
24     from esys.escript import *
25     from esys.pycad import *
26     from esys.pycad.gmsh import Design
27     from esys.escript.unitsSI import *
28     from math import tan
29 sshaw 5288
30     try:
31     # This imports the rectangle domain function
32     from esys.finley import MakeDomain
33     HAVE_FINLEY = True
34     except ImportError:
35     print("Finley module not available")
36     HAVE_FINLEY = False
37 gross 3439 #the folder to put our outputs in, leave blank "" for script path
38     save_path= os.path.join("data","example10")
39 ahallam 3135 #
40     # input data:
41     #
42     width=5*km
43     length=5*km
44     depth=2*km
45     fault_w_front=200*m
46     fault_w_back=200*m
47     fault_mid_front=2*km
48     fault_mid_back=4*km
49 ahallam 3195 fault_dip_front=30*DEG
50     fault_dip_back=30*DEG
51 ahallam 3135
52 ahallam 3195 layers_left_at_front=[ [ 'limestone' , 1100*m ] , [ 'xx' , 150*m ] ,['shale',300*m], [ 'limestone' ] ]
53     layers_right_at_front=[ [ 'limestone' , 500*m ], [ 'xx' , 150*m ],['shale',300*m], ['limestone' ] ]
54 ahallam 3135 slop=2*DEG
55     #
56     # ====================================================================================
57     #
58     def setOffsetsFromThickness(layers):
59     l=0
60     offsets=[]
61     for n in layers:
62     if len(n)<2:
63     d=depth-l
64     if depth-l <=0:
65     raise "Layer structure is too thick. Increase depth > %e"%l
66     else:
67     d=n[1]
68     l+=d
69     offsets.append([ n[0], -(l-d), d, -l ])
70     return offsets
71     def setVerticalPositionsFromDip(layers, dip):
72     offsets=[]
73     s=1/tan(90*DEG-dip)
74     for n in layers:
75     name=n[0]
76     top=n[1]
77     d=n[2]
78     bot=n[3]
79     offsets.append([ name, top, s*top, d, bot, s*bot ] )
80     return offsets
81    
82     def setBackLayers(layers, slop, dip):
83     if layers[0][3]-tan(slop)*length <=0:
84 jfenwick 3892 raise ValueError("negative thickness in %s at back"%(layers[0][0],))
85 ahallam 3135 out=[[ layers[0][0], layers[0][3]-tan(slop)*length ]]
86     for i in range(1,len(layers)-1): out.append([layers[i][0], layers[i][3] ])
87     out.append([ layers[-1][0] ])
88     out=setOffsetsFromThickness(out)
89     return setVerticalPositionsFromDip(out, dip)
90    
91     def getCutLine(p0,layers,offset=True):
92     out=[ ]
93     p=p0
94     for n in layers:
95     if offset:
96     p1,p = p, Point(p0.getCoordinates()[0]+n[5],p0.getCoordinates()[1],p0.getCoordinates()[2]+n[4])
97     else:
98     p1,p = p, Point(p0.getCoordinates()[0],p0.getCoordinates()[1],p0.getCoordinates()[2]+n[4])
99     out.append(Line(p1,p))
100     return out
101    
102     def getBackLine(front, back):
103     out = [ Line(front[0].getStartPoint(), back[0].getStartPoint()) ]
104     for i in range(len(front)):
105     out.append(Line(front[i].getEndPoint(), back[i].getEndPoint()))
106     return out
107    
108     def getCrossLine(left, right):
109     return getBackLine(left,right)
110    
111     def addVolume(front_left, back_left, front_right, back_right, PS, FF, map, filter_left=False):
112    
113     front_to_back_left = getBackLine(front_left , back_left)
114     front_to_back_right = getBackLine(front_right , back_right)
115     front_left_to_right= getCrossLine(front_left, front_right)
116     back_left_to_right= getCrossLine(back_left, back_right)
117    
118     #if filter_left:
119     # out1=front_to_back_left[0]
120     # out2=front_to_back_left[-1]
121     #else:
122     out1=front_to_back_right[0]
123     out2=front_to_back_right[-1]
124    
125     topface=PlaneSurface(CurveLoop(front_left_to_right[0], front_to_back_right[0], -back_left_to_right[0], -front_to_back_left[0]))
126     for i in range(len(map)):
127     name=map[i][0]
128     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]))
129     face_front=PlaneSurface(CurveLoop(front_left_to_right[i+1], -front_right[i], -front_left_to_right[i], front_left[i]))
130     face_back=PlaneSurface(CurveLoop(back_left_to_right[i+1], -back_right[i], -back_left_to_right[i], back_left[i]))
131     face_left=PlaneSurface(CurveLoop(front_left[i], front_to_back_left[i+1], -back_left[i], -front_to_back_left[i]))
132     face_right=PlaneSurface(CurveLoop(front_right[i], front_to_back_right[i+1], -back_right[i], -front_to_back_right[i]))
133    
134     v=Volume(SurfaceLoop(topface,-face2, face_front, -face_back, -face_left, face_right))
135 jfenwick 3892 print(v)
136 ahallam 3135 topface=face2
137     if filter_left:
138     FF.append(face_right)
139     else:
140     FF.append(-face_right)
141 jfenwick 3892 if name not in PS: PS[name]=PropertySet(name)
142 ahallam 3135 PS[name].addItem(v)
143     return out1, out2, PS, FF
144    
145    
146    
147 sshaw 5288 if HAVE_FINLEY:
148     layers_left_at_edge_at_front=setOffsetsFromThickness(layers_left_at_front)
149     layers_right_at_edge_at_front=setOffsetsFromThickness(layers_right_at_front)
150 ahallam 3135
151 sshaw 5288 layers_left_at_front=setVerticalPositionsFromDip(layers_left_at_edge_at_front, fault_dip_front)
152     layers_right_at_front=setVerticalPositionsFromDip(layers_right_at_edge_at_front, fault_dip_front)
153 ahallam 3135
154    
155 sshaw 5288 layers_left_at_back=setBackLayers(layers_left_at_front, slop, fault_dip_back)
156     layers_right_at_back=setBackLayers(layers_right_at_front, slop, fault_dip_back)
157 ahallam 3135
158    
159 sshaw 5288 left_front_edge=getCutLine(Point(0.,0.,0.), layers_left_at_front, offset=False)
160     left_front_fault=getCutLine(Point(fault_mid_front-fault_w_front/2,0.,0.), layers_left_at_front, offset=True)
161 ahallam 3135
162 sshaw 5288 right_front_fault=getCutLine(Point(fault_mid_front+fault_w_front/2,0.,0.), layers_right_at_front, offset=True)
163     right_front_edge=getCutLine(Point(width,0.,0.), layers_right_at_front, offset=False)
164 ahallam 3135
165 sshaw 5288 left_back_edge=getCutLine(Point(0.,length,0.), layers_left_at_back, offset=False)
166     left_back_fault=getCutLine(Point(fault_mid_back-fault_w_back/2,length,0.), layers_left_at_back, offset=True)
167 ahallam 3135
168 sshaw 5288 right_back_fault=getCutLine(Point(fault_mid_back+fault_w_back/2,length,0.), layers_right_at_back, offset=True)
169     right_back_edge=getCutLine(Point(width,length,0.), layers_right_at_back, offset=False)
170 ahallam 3135
171 sshaw 5288 PS={}
172     FF=[]
173     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)
174     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)
175 ahallam 3135
176 sshaw 5288 fault_line_top_front=Line(front_to_back_left_top.getStartPoint(), front_to_back_right_top.getStartPoint())
177     fault_line_bot_front=Line(front_to_back_left_bot.getStartPoint(), front_to_back_right_bot.getStartPoint())
178     fault_line_top_back=Line(front_to_back_left_top.getEndPoint(), front_to_back_right_top.getEndPoint())
179     fault_line_bot_back=Line(front_to_back_left_bot.getEndPoint(), front_to_back_right_bot.getEndPoint())
180 ahallam 3135
181 sshaw 5288 FF.append(PlaneSurface(CurveLoop(front_to_back_left_top,fault_line_top_back,-front_to_back_right_top,-fault_line_top_front)))
182     FF.append(-PlaneSurface(CurveLoop(front_to_back_left_bot,fault_line_bot_back,-front_to_back_right_bot,-fault_line_bot_front)))
183 ahallam 3135
184 sshaw 5288 FF.append(PlaneSurface(CurveLoop(*tuple([ -fault_line_top_front,fault_line_bot_front ]+left_front_fault+[ -l for l in right_front_fault ]))))
185     FF.append(-PlaneSurface(CurveLoop(*tuple([ -fault_line_top_back,fault_line_bot_back ]+left_back_fault+[ -l for l in right_back_fault ]))))
186 ahallam 3135
187    
188 sshaw 5288 # war 120
189     des=Design(dim=3, order=1, element_size = 400*m, keep_files=True)
190     des.addItems(*tuple(PS.values()))
191     des.addItems(PropertySet("fault",Volume(SurfaceLoop( *tuple(FF)))))
192     des.setMeshFileName(os.path.join(save_path,"fault.msh"))
193     dom=MakeDomain(des)
194     dom.write(os.path.join(save_path,"fault.fly"))
195 ahallam 3135

  ViewVC Help
Powered by ViewVC 1.1.26