| 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 |
|
Parent Directory
| 
Revision Log