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