/[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 3892 - (hide annotations)
Tue Apr 10 08:57:23 2012 UTC (7 years ago) by jfenwick
File MIME type: text/x-python
File size: 7850 byte(s)
Merged changes across from the attempt2 branch.
This version builds and passes python2 tests.
It also passes most python3 tests.



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

  ViewVC Help
Powered by ViewVC 1.1.26