| 1 |
# $Id$ |
| 2 |
|
| 3 |
""" |
| 4 |
|
| 5 |
generates finley mesh simple vertical fault |
| 6 |
|
| 7 |
THIS CODE CREATES RICH CONTACT ELEMENTS AND RICH FACE ELEMENTS |
| 8 |
with fix for contact elements at FAULT ENDS |
| 9 |
|
| 10 |
|
| 11 |
@var __author__: name of author |
| 12 |
@var __copyright__: copyrights |
| 13 |
@var __license__: licence agreement |
| 14 |
@var __url__: url entry point on documentation |
| 15 |
@var __version__: version |
| 16 |
@var __date__: date of the version |
| 17 |
""" |
| 18 |
|
| 19 |
__author__="Louise Kettle" |
| 20 |
__copyright__=""" Copyright (c) 2006 by ACcESS MNRF |
| 21 |
http://www.access.edu.au |
| 22 |
Primary Business: Queensland, Australia""" |
| 23 |
__license__="""Licensed under the Open Software License version 3.0 |
| 24 |
http://www.opensource.org/licenses/osl-3.0.php""" |
| 25 |
__url__="http://www.iservo.edu.au/esys" |
| 26 |
__version__="$Revision$" |
| 27 |
__date__="$Date$" |
| 28 |
|
| 29 |
from esys.escript import * |
| 30 |
from numarray import zeros,Float,array,size |
| 31 |
|
| 32 |
#... generate domain ... |
| 33 |
ne = 18 |
| 34 |
width = 100000. |
| 35 |
height = 30000. |
| 36 |
fstart=numarray.array([width/2.,7.*width/16.,3.*height/8.]) |
| 37 |
fend=numarray.array([width/2.,9.*width/16.,5.*height/8.]) |
| 38 |
|
| 39 |
def faultL(l0,l1, l2,ne0, ne1, ne2,contact=False,xstart=zeros(3),xend=zeros(3)): |
| 40 |
# if contact = true then there is a fault surface |
| 41 |
# xstart is the co-ordinates of start of fault |
| 42 |
# xend is the co-ordinates at end of fault |
| 43 |
meshfaultL=open('meshfault3D.fly','w') |
| 44 |
|
| 45 |
FaultError1="ERROR: fault defined on or too close to an outer surface" |
| 46 |
FaultError2="ERROR: the mesh is too coarse for fault" |
| 47 |
|
| 48 |
N0=ne0+1 |
| 49 |
N1=ne1+1 |
| 50 |
N2=ne2+1 |
| 51 |
if (N0<=N1 and N0<=N2): |
| 52 |
if (N1 <= N2): |
| 53 |
M0=1 |
| 54 |
M1=N0 |
| 55 |
M2=N0*N1 |
| 56 |
M0i=1 |
| 57 |
M1i=ne0 |
| 58 |
M2i=ne0*ne1 |
| 59 |
else: |
| 60 |
M0=1 |
| 61 |
M2=N0 |
| 62 |
M1=N0*N2 |
| 63 |
M0i=1 |
| 64 |
M2i=ne0 |
| 65 |
M1i=ne0*ne2 |
| 66 |
elif (N1<=N2 and N1<=N0): |
| 67 |
if (N2 <= N0): |
| 68 |
M1=1 |
| 69 |
M2=N1 |
| 70 |
M0=N2*N1 |
| 71 |
M1i=1 |
| 72 |
M2i=ne1 |
| 73 |
M0i=ne2*ne1 |
| 74 |
else: |
| 75 |
M1=1 |
| 76 |
M0=N1 |
| 77 |
M2=N1*N0 |
| 78 |
M1i=1 |
| 79 |
M0i=ne1 |
| 80 |
M2i=ne0*ne1 |
| 81 |
else: |
| 82 |
if (N0 <= N1): |
| 83 |
M2=1 |
| 84 |
M0=N2 |
| 85 |
M1=N2*N0 |
| 86 |
M2i=1 |
| 87 |
M0i=ne2 |
| 88 |
M1i=ne0*ne2 |
| 89 |
else: |
| 90 |
M2=1 |
| 91 |
M1=N2 |
| 92 |
M0=N1*N2 |
| 93 |
M2i=1 |
| 94 |
M1i=ne2 |
| 95 |
M0i=ne2*ne1 |
| 96 |
|
| 97 |
dim=3 |
| 98 |
Element_numNodes=8 |
| 99 |
Element_Num=ne0*ne1*ne2 |
| 100 |
if contact==False: |
| 101 |
numNodes=N0*N1*N2 |
| 102 |
|
| 103 |
else: |
| 104 |
# define double (contact element) nodes on interior of fault |
| 105 |
i0start=round(xstart[0]*ne0/l0) |
| 106 |
i1start=round(xstart[1]*ne1/l1) |
| 107 |
i2start=round(xstart[2]*ne2/l2) |
| 108 |
i0end=round(xend[0]*ne0/l0) |
| 109 |
i1end=round(xend[1]*ne1/l1) |
| 110 |
i2end=round(xend[2]*ne2/l2) |
| 111 |
n0double=int(i0end)-int(i0start) |
| 112 |
n1double=int(i1end)-int(i1start) |
| 113 |
n2double=int(i2end)-int(i2start) |
| 114 |
print "fstart = ",[i0start*l0/ne0, i1start*l1/ne1, i2start*l2/ne2] |
| 115 |
print "fend = ", [i0start*l0/ne0, i1start*l1/ne1 + n1double*l1/ne1,i2start*l2/ne2 + n2double*l2/ne2 -height] |
| 116 |
if (i0start == 0) or (i1start==0) or (i2start==0): |
| 117 |
raise FaultError1 |
| 118 |
|
| 119 |
if (i0end == ne0) or (i1end==ne1) or (i2end==ne2): |
| 120 |
raise FaultError1 |
| 121 |
|
| 122 |
if n0double==0: |
| 123 |
numNodes=N0*N1*N2+(n1double-1)*(n2double-1) |
| 124 |
|
| 125 |
elif n1double==0: |
| 126 |
numNodes=N0*N1*N2+(n0double-1)*(n2double-1) |
| 127 |
|
| 128 |
elif n2double==0: |
| 129 |
numNodes=N0*N1*N2+(n0double-1)*(n1double-1) |
| 130 |
|
| 131 |
# define nodes for normal elements |
| 132 |
# there are N0*N1*N2 normal nodes |
| 133 |
|
| 134 |
Node=zeros([3,numNodes],Float) |
| 135 |
Node_ref=zeros(numNodes,Float) |
| 136 |
Node_DOF=zeros(numNodes,Float) |
| 137 |
Node_tag=zeros(numNodes,Float) |
| 138 |
|
| 139 |
meshfaultL.write("KettleFault\n") |
| 140 |
#print 'Nodes' |
| 141 |
meshfaultL.write("%dD-nodes %d\n"%(dim,numNodes)) |
| 142 |
|
| 143 |
for i2 in range(N2): |
| 144 |
for i1 in range (N1): |
| 145 |
for i0 in range(N0): |
| 146 |
k=M0*i0+M1*i1+M2*i2; |
| 147 |
Node_ref[k]=i0 + N0*i1 + N0*N1*i2 |
| 148 |
# no periodic boundary conditions |
| 149 |
Node_DOF[k]=Node_ref[k] |
| 150 |
Node_tag[k]=0 |
| 151 |
Node[0][k]=(i0)*l0/(N0-1) |
| 152 |
Node[1][k]=(i1)*l1/(N1-1) |
| 153 |
Node[2][k]=(i2)*l2/(N2-1)-height |
| 154 |
|
| 155 |
# define double nodes on fault (will have same coordinates as some of nodes already defined) |
| 156 |
# only get double nodes on INTERIOR of fault |
| 157 |
|
| 158 |
if contact==True: |
| 159 |
Fault_NE=N0*N1*N2 |
| 160 |
if n0double==0: |
| 161 |
if(n1double<=n2double): |
| 162 |
M1f=1 |
| 163 |
M2f=n1double-1 |
| 164 |
else: |
| 165 |
M2f=1 |
| 166 |
M1f=n2double-1 |
| 167 |
|
| 168 |
for i2 in range(n2double-1): |
| 169 |
for i1 in range(n1double-1): |
| 170 |
k=Fault_NE+M1f*i1+M2f*i2 |
| 171 |
Node_ref[k]=Fault_NE + i1 + (n1double-1)*i2 |
| 172 |
Node_DOF[k]=Node_ref[k] |
| 173 |
Node_tag[k]=1 |
| 174 |
Node[0][k]=i0start*l0/ne0 |
| 175 |
Node[1][k]=i1start*l1/ne1 + (i1+1)*l1/ne1 |
| 176 |
Node[2][k]=i2start*l2/ne2 + (i2+1)*l2/ne2 -height |
| 177 |
|
| 178 |
# elif n1double==0: |
| 179 |
# elif n2double==0: |
| 180 |
|
| 181 |
# write nodes to file |
| 182 |
for i in range(numNodes): |
| 183 |
meshfaultL.write("%d %d %d"%(Node_ref[i],Node_DOF[i],Node_tag[i])) |
| 184 |
for j in range(dim): |
| 185 |
meshfaultL.write(" %lf"%Node[j][i]) |
| 186 |
meshfaultL.write("\n") |
| 187 |
|
| 188 |
|
| 189 |
|
| 190 |
|
| 191 |
# defining interior elements |
| 192 |
# there are ne0*ne1*ne2 interior elements |
| 193 |
|
| 194 |
Element_Nodes=zeros([8,ne0*ne1*ne2],Float) |
| 195 |
Element_ref=zeros(ne0*ne1*ne2,Float) |
| 196 |
Element_tag=zeros(ne0*ne1*ne2,Float) |
| 197 |
|
| 198 |
#print 'Interior elements' |
| 199 |
|
| 200 |
for i2 in range(ne2): |
| 201 |
for i1 in range (ne1): |
| 202 |
for i0 in range(ne0): |
| 203 |
k=i0 + ne0*i1 + ne0*ne1*i2; |
| 204 |
# define corner node (node0) |
| 205 |
node0=i0 + N0*i1 + N0*N1*i2; |
| 206 |
Element_ref[k]=k |
| 207 |
Element_tag[k]=0 |
| 208 |
# for hex8 the interior elements are specified by 8 nodes |
| 209 |
Element_Nodes[0][k]=node0; |
| 210 |
Element_Nodes[1][k]=node0+1; |
| 211 |
Element_Nodes[2][k]=node0+N0+1; |
| 212 |
Element_Nodes[3][k]=node0+N0; |
| 213 |
Element_Nodes[4][k]=node0+N0*N1; |
| 214 |
Element_Nodes[5][k]=node0+N0*N1+1; |
| 215 |
Element_Nodes[6][k]=node0+N0*N1+N0+1; |
| 216 |
Element_Nodes[7][k]=node0+N0*N1+N0; |
| 217 |
|
| 218 |
if contact==True: |
| 219 |
if n0double==0: |
| 220 |
x0s= i0start*l0/ne0 |
| 221 |
x1s= i1start*l1/ne1 |
| 222 |
x2s= i2start*l2/ne2 |
| 223 |
x0e= i0end*l0/ne0 |
| 224 |
x1e= i1end*l1/ne1 |
| 225 |
x2e= i2end*l2/ne2 |
| 226 |
#print x0s,x1s,x2s,x0e,x1e,x2e |
| 227 |
|
| 228 |
if (n1double==1) or (n2double==1): |
| 229 |
raise FaultError2 |
| 230 |
for i2 in range(n2double): |
| 231 |
for i1 in range(n1double): |
| 232 |
# here the coordinates of kfault and kold are the same |
| 233 |
# Ref for fault node (only on interior nodes of fault): |
| 234 |
if (i1>0) and (i2>0): |
| 235 |
kfault=Fault_NE+M1f*(i1-1.)+M2f*(i2-1.) |
| 236 |
#print kfault, Node[0][kfault],Node[1][kfault],Node[2][kfault] |
| 237 |
else: |
| 238 |
kfault=0. |
| 239 |
# determine bottom corner node of each element |
| 240 |
|
| 241 |
# Ref for normal interior node: |
| 242 |
kold=int(M0*i0start+M1*(i1start + i1) + M2*(i2start+i2)) |
| 243 |
#print kold, Node[0][kold],Node[1][kold],Node[2][kold] |
| 244 |
# Ref for interior element: |
| 245 |
kint=int(M0i*i0start + M1i*(i1start+i1) + M2i*(i2start+i2)) |
| 246 |
#print kint, Element_Nodes[0][kint] |
| 247 |
|
| 248 |
x0= (i0start)*l0/ne0 |
| 249 |
x1= (i1start+i1)*l1/ne1 |
| 250 |
x2= (i2start+i2)*l2/ne2 |
| 251 |
|
| 252 |
# for x0 > xstart we need to overwrite old Nodes in interior element references |
| 253 |
# with fault nodes: |
| 254 |
|
| 255 |
# for the interior elements with x1<x1s and x2<x2s the only nodes need changing |
| 256 |
# are on the fault: |
| 257 |
if (i1==0) and (i2==0): |
| 258 |
# nearest fault node: |
| 259 |
kfaultref=int(Fault_NE+M1f*i1+M2f*i2) |
| 260 |
elif (i1==0): |
| 261 |
# nearest fault node |
| 262 |
kfaultref=int(Fault_NE+i1*M1f+(i2-1.)*(M2f)) |
| 263 |
elif (i2==0): |
| 264 |
# nearest fault node |
| 265 |
kfaultref=int(Fault_NE+(i1-1.)*M1f + i2*M2f) |
| 266 |
else: |
| 267 |
# looking at element with fault node on bottom corner |
| 268 |
kfaultref=int(kfault) |
| 269 |
|
| 270 |
#print x0,x1,x2 |
| 271 |
#print kold, Node[0][kold],Node[1][kold],Node[2][kold] |
| 272 |
#print kfaultref, Node[0][kfaultref],Node[1][kfaultref],Node[2][kfaultref] |
| 273 |
|
| 274 |
# overwrite 4 outer corner elements of fault (only one node changed) |
| 275 |
if (i1==0 and i2==0): |
| 276 |
Element_Nodes[7][kint]=kfaultref |
| 277 |
|
| 278 |
elif (i1==0 and i2==n2double-1): |
| 279 |
Element_Nodes[3][kint]=kfaultref |
| 280 |
|
| 281 |
elif (i1==n1double-1 and i2==0): |
| 282 |
Element_Nodes[4][kint]=kfaultref |
| 283 |
|
| 284 |
elif (i1==n1double-1 and i2==n2double-1): |
| 285 |
Element_Nodes[0][kint]=kfaultref |
| 286 |
|
| 287 |
# overwrite 4 sides of fault (only 2 nodes changed) |
| 288 |
elif (i1==0): |
| 289 |
Element_Nodes[3][kint]=kfaultref |
| 290 |
kfaultref1=int(kfaultref+M2f) |
| 291 |
Element_Nodes[7][kint]=kfaultref1 |
| 292 |
|
| 293 |
elif (i1==n1double-1): |
| 294 |
Element_Nodes[0][kint]=kfaultref |
| 295 |
kfaultref1=kfaultref+M2f |
| 296 |
Element_Nodes[4][kint]=kfaultref1 |
| 297 |
|
| 298 |
elif (i2==0): |
| 299 |
Element_Nodes[4][kint]=kfaultref |
| 300 |
kfaultref1=kfaultref+M1f |
| 301 |
Element_Nodes[7][kint]=kfaultref1 |
| 302 |
|
| 303 |
elif (i2==n2double-1): |
| 304 |
#print i1 |
| 305 |
#nodecheck=int(Element_Nodes[0][kint] ) |
| 306 |
#print nodecheck, Node[0][nodecheck],Node[1][nodecheck],Node[2][nodecheck] |
| 307 |
#print kfaultref, Node[0][kfaultref],Node[1][kfaultref],Node[2][kfaultref] |
| 308 |
|
| 309 |
Element_Nodes[0][kint]=kfaultref |
| 310 |
kfaultref1=kfaultref+M1f |
| 311 |
|
| 312 |
#nodecheck=int(Element_Nodes[3][kint] ) |
| 313 |
#print nodecheck, Node[0][nodecheck],Node[1][nodecheck],Node[2][nodecheck] |
| 314 |
#print kfaultref1, Node[0][kfaultref1],Node[1][kfaultref1],Node[2][kfaultref1] |
| 315 |
|
| 316 |
Element_Nodes[3][kint]=kfaultref1 |
| 317 |
|
| 318 |
|
| 319 |
# overwrite interior fault elements (4 nodes changed) |
| 320 |
else: |
| 321 |
#print i1,i2 |
| 322 |
#nodecheck=int(Element_Nodes[0][kint] ) |
| 323 |
#print nodecheck, Node[0][nodecheck],Node[1][nodecheck],Node[2][nodecheck] |
| 324 |
#print kfaultref, Node[0][kfaultref],Node[1][kfaultref],Node[2][kfaultref] |
| 325 |
|
| 326 |
Element_Nodes[0][kint]=kfaultref |
| 327 |
#if (x1<x1e and x2<x2e): |
| 328 |
kfaultref1=kfaultref+M1f |
| 329 |
kfaultref2=kfaultref+M2f |
| 330 |
kfaultref3=kfaultref+M1f+M2f |
| 331 |
|
| 332 |
nodecheck=int(Element_Nodes[3][kint] ) |
| 333 |
#print nodecheck, Node[0][nodecheck],Node[1][nodecheck],Node[2][nodecheck] |
| 334 |
#print kfaultref1, Node[0][kfaultref1],Node[1][kfaultref1],Node[2][kfaultref1] |
| 335 |
|
| 336 |
nodecheck=int(Element_Nodes[4][kint] ) |
| 337 |
#print nodecheck, Node[0][nodecheck],Node[1][nodecheck],Node[2][nodecheck] |
| 338 |
#print kfaultref2, Node[0][kfaultref2],Node[1][kfaultref2],Node[2][kfaultref2] |
| 339 |
|
| 340 |
nodecheck=int(Element_Nodes[7][kint] ) |
| 341 |
#print nodecheck, Node[0][nodecheck],Node[1][nodecheck],Node[2][nodecheck] |
| 342 |
#print kfaultref3, Node[0][kfaultref3],Node[1][kfaultref3],Node[2][kfaultref3] |
| 343 |
|
| 344 |
|
| 345 |
Element_Nodes[3][kint]=kfaultref1 |
| 346 |
Element_Nodes[4][kint]=kfaultref2 |
| 347 |
Element_Nodes[7][kint]=kfaultref3 |
| 348 |
#elif x1<x1e: |
| 349 |
# kfaultref1=kfaultref+M1f |
| 350 |
# Element_Nodes[3][kint]=kfaultref1 |
| 351 |
#elif x2<x2e: |
| 352 |
# kfaultref2=kfaultref+M2f |
| 353 |
# Element_Nodes[4][kint]=kfaultref2 |
| 354 |
|
| 355 |
|
| 356 |
# print kint, kfaultref |
| 357 |
|
| 358 |
# write interior elements to file |
| 359 |
Element_Type='Hex8' |
| 360 |
meshfaultL.write("%s %d\n"%(Element_Type,Element_Num)) |
| 361 |
|
| 362 |
for i in range(Element_Num): |
| 363 |
meshfaultL.write("%d %d"%(Element_ref[i],Element_tag[i])) |
| 364 |
for j in range(Element_numNodes): |
| 365 |
meshfaultL.write(" %d"%Element_Nodes[j][i]) |
| 366 |
meshfaultL.write("\n") |
| 367 |
|
| 368 |
# face elements |
| 369 |
FaceElement_Type='Hex8Face' |
| 370 |
FaceElement_Num= 2*(ne0*ne1 + ne0*ne2 + ne1*ne2) |
| 371 |
FaceElement_numNodes=8 |
| 372 |
|
| 373 |
meshfaultL.write("%s %d\n"%(FaceElement_Type,FaceElement_Num)) |
| 374 |
|
| 375 |
FaceElement_Nodes=zeros([FaceElement_numNodes,FaceElement_Num],Float) |
| 376 |
FaceElement_ref=zeros(FaceElement_Num,Float) |
| 377 |
FaceElement_tag=zeros(FaceElement_Num,Float) |
| 378 |
|
| 379 |
kcount=0 |
| 380 |
|
| 381 |
# defining face elements on x2=0 face |
| 382 |
for i1 in range (ne1): |
| 383 |
for i0 in range(ne0): |
| 384 |
i2=0 |
| 385 |
k=i0 + ne0*i1 + ne0*ne1*i2; |
| 386 |
# define corner node (node0) |
| 387 |
node0=i0 + N0*i1 + N0*N1*i2; |
| 388 |
FaceElement_ref[kcount]=kcount |
| 389 |
FaceElement_tag[kcount]=3 |
| 390 |
# for hex8face the face elements are specified by 8 nodes |
| 391 |
FaceElement_Nodes[0][kcount]=node0; |
| 392 |
FaceElement_Nodes[1][kcount]=node0+1; |
| 393 |
FaceElement_Nodes[2][kcount]=node0+N0+1; |
| 394 |
FaceElement_Nodes[3][kcount]=node0+N0; |
| 395 |
FaceElement_Nodes[4][kcount]=node0+N0*N1; |
| 396 |
FaceElement_Nodes[5][kcount]=node0+N0*N1+1; |
| 397 |
FaceElement_Nodes[6][kcount]=node0+N0*N1+N0+1; |
| 398 |
FaceElement_Nodes[7][kcount]=node0+N0*N1+N0; |
| 399 |
kcount+=1 |
| 400 |
|
| 401 |
# defining face elements on x2=L face |
| 402 |
for i1 in range (ne1): |
| 403 |
for i0 in range(ne0): |
| 404 |
i2=ne2-1 |
| 405 |
k=i0 + ne0*i1 + ne0*ne1*i2; |
| 406 |
# define corner node (node0) |
| 407 |
node0=i0 + N0*i1 + N0*N1*i2; |
| 408 |
FaceElement_ref[kcount]=kcount |
| 409 |
FaceElement_tag[kcount]=3 |
| 410 |
# for hex8face the face elements are specified by 8 nodes |
| 411 |
FaceElement_Nodes[0][kcount]=node0+N0*N1; |
| 412 |
FaceElement_Nodes[1][kcount]=node0+N0*N1+1; |
| 413 |
FaceElement_Nodes[2][kcount]=node0+N0*N1+N0+1; |
| 414 |
FaceElement_Nodes[3][kcount]=node0+N0*N1+N0; |
| 415 |
FaceElement_Nodes[4][kcount]=node0; |
| 416 |
FaceElement_Nodes[5][kcount]=node0+1; |
| 417 |
FaceElement_Nodes[6][kcount]=node0+N0+1; |
| 418 |
FaceElement_Nodes[7][kcount]=node0+N0; |
| 419 |
kcount+=1 |
| 420 |
|
| 421 |
# defining face elements on x1=0 face |
| 422 |
for i2 in range (ne2): |
| 423 |
for i0 in range(ne0): |
| 424 |
i1=0 |
| 425 |
k=i0 + ne0*i1 + ne0*ne1*i2; |
| 426 |
# define corner node (node0) |
| 427 |
node0=i0 + N0*i1 + N0*N1*i2; |
| 428 |
FaceElement_ref[kcount]=kcount |
| 429 |
FaceElement_tag[kcount]=3 |
| 430 |
# for hex8face the face elements are specified by 8 nodes |
| 431 |
FaceElement_Nodes[0][kcount]=node0; |
| 432 |
FaceElement_Nodes[1][kcount]=node0+N0*N1; |
| 433 |
FaceElement_Nodes[2][kcount]=node0+N0*N1+1; |
| 434 |
FaceElement_Nodes[3][kcount]=node0+1; |
| 435 |
FaceElement_Nodes[4][kcount]=node0+N0; |
| 436 |
FaceElement_Nodes[5][kcount]=node0+N0*N1+N0; |
| 437 |
FaceElement_Nodes[6][kcount]=node0+N0*N1+N0+1; |
| 438 |
FaceElement_Nodes[7][kcount]=node0+N0+1; |
| 439 |
kcount+=1 |
| 440 |
|
| 441 |
# defining face elements on x1=L face |
| 442 |
for i2 in range (ne2): |
| 443 |
for i0 in range(ne0): |
| 444 |
i1=ne1-1 |
| 445 |
k=i0 + ne0*i1 + ne0*ne1*i2; |
| 446 |
# define corner node (node0) |
| 447 |
node0=i0 + N0*i1 + N0*N1*i2; |
| 448 |
FaceElement_ref[kcount]=kcount |
| 449 |
FaceElement_tag[kcount]=3 |
| 450 |
# for hex8face the face elements are specified by 8 nodes |
| 451 |
FaceElement_Nodes[0][kcount]=node0+N0; |
| 452 |
FaceElement_Nodes[1][kcount]=node0+N0*N1+N0; |
| 453 |
FaceElement_Nodes[2][kcount]=node0+N0*N1+N0+1; |
| 454 |
FaceElement_Nodes[3][kcount]=node0+N0+1; |
| 455 |
FaceElement_Nodes[4][kcount]=node0; |
| 456 |
FaceElement_Nodes[5][kcount]=node0+N0*N1; |
| 457 |
FaceElement_Nodes[6][kcount]=node0+N0*N1+1; |
| 458 |
FaceElement_Nodes[7][kcount]=node0+1; |
| 459 |
kcount+=1 |
| 460 |
|
| 461 |
# defining face elements on x0=0 face |
| 462 |
for i2 in range (ne2): |
| 463 |
for i1 in range(ne1): |
| 464 |
i0=0 |
| 465 |
k=i0 + ne0*i1 + ne0*ne1*i2; |
| 466 |
# define corner node (node0) |
| 467 |
node0=i0 + N0*i1 + N0*N1*i2; |
| 468 |
FaceElement_ref[kcount]=kcount |
| 469 |
FaceElement_tag[kcount]=3 |
| 470 |
# for hex8face the face elements are specified by 8 nodes |
| 471 |
FaceElement_Nodes[0][kcount]=node0; |
| 472 |
FaceElement_Nodes[1][kcount]=node0+N0; |
| 473 |
FaceElement_Nodes[2][kcount]=node0+N0*N1+N0; |
| 474 |
FaceElement_Nodes[3][kcount]=node0+N0*N1; |
| 475 |
FaceElement_Nodes[4][kcount]=node0+1; |
| 476 |
FaceElement_Nodes[5][kcount]=node0+N0+1; |
| 477 |
FaceElement_Nodes[6][kcount]=node0+N0*N1+N0+1; |
| 478 |
FaceElement_Nodes[7][kcount]=node0+N0*N1+1; |
| 479 |
kcount+=1 |
| 480 |
|
| 481 |
# defining face elements on x0=L face |
| 482 |
for i2 in range (ne2): |
| 483 |
for i1 in range(ne1): |
| 484 |
i0=ne1-1 |
| 485 |
k=i0 + ne0*i1 + ne0*ne1*i2; |
| 486 |
# define corner node (node0) |
| 487 |
node0=i0 + N0*i1 + N0*N1*i2; |
| 488 |
FaceElement_ref[kcount]=kcount |
| 489 |
FaceElement_tag[kcount]=3 |
| 490 |
# for hex8face the face elements are specified by 8 nodes |
| 491 |
FaceElement_Nodes[0][kcount]=node0+1; |
| 492 |
FaceElement_Nodes[1][kcount]=node0+N0+1; |
| 493 |
FaceElement_Nodes[2][kcount]=node0+N0*N1+N0+1; |
| 494 |
FaceElement_Nodes[3][kcount]=node0+N0*N1+1; |
| 495 |
FaceElement_Nodes[4][kcount]=node0; |
| 496 |
FaceElement_Nodes[5][kcount]=node0+N0; |
| 497 |
FaceElement_Nodes[6][kcount]=node0+N0*N1+N0; |
| 498 |
FaceElement_Nodes[7][kcount]=node0+N0*N1; |
| 499 |
kcount+=1 |
| 500 |
|
| 501 |
|
| 502 |
|
| 503 |
|
| 504 |
for i in range(FaceElement_Num): |
| 505 |
meshfaultL.write("%d %d"%(FaceElement_ref[i],FaceElement_tag[i])) |
| 506 |
for j in range(FaceElement_numNodes): |
| 507 |
meshfaultL.write(" %d"%FaceElement_Nodes[j][i]) |
| 508 |
meshfaultL.write("\n") |
| 509 |
|
| 510 |
|
| 511 |
# contact elements |
| 512 |
ContactElement_Type='Hex8Face_Contact' |
| 513 |
ContactElement_Num=0 |
| 514 |
ContactElement_numNodes=16 |
| 515 |
# print contact elements on fault |
| 516 |
if contact==True: |
| 517 |
if n0double==0: |
| 518 |
ContactElement_Num=(n1double)*(n2double) |
| 519 |
ContactElement_Nodes=zeros([ContactElement_numNodes,ContactElement_Num],Float) |
| 520 |
ContactElement_ref=zeros(ContactElement_Num,Float) |
| 521 |
ContactElement_tag=zeros(ContactElement_Num,Float) |
| 522 |
#print ContactElement_Num |
| 523 |
|
| 524 |
for i2 in range(n2double): |
| 525 |
for i1 in range(n1double): |
| 526 |
k=i1+(n1double)*i2 |
| 527 |
#print k |
| 528 |
# define reference for interior elements with x0<=x0s |
| 529 |
# here the nodes are the old interior nodes |
| 530 |
kintold=int(M0i*(i0start-1) + M1i*(i1start+i1) + M2i*(i2start+i2)) |
| 531 |
|
| 532 |
# define reference for interior elements with x0>x0s |
| 533 |
# here the double nodes are the fault nodes |
| 534 |
|
| 535 |
kintfault=int(M0i*i0start + M1i*(i1start+i1) + M2i*(i2start+i2)) |
| 536 |
|
| 537 |
#nodecheck=int(Element_Nodes[1][kintold] ) |
| 538 |
#print nodecheck, Node[0][nodecheck],Node[1][nodecheck],Node[2][nodecheck] |
| 539 |
#nodecheck=int(Element_Nodes[0][kintfault] ) |
| 540 |
#print nodecheck, Node[0][nodecheck],Node[1][nodecheck],Node[2][nodecheck] |
| 541 |
|
| 542 |
#nodecheck=int(Element_Nodes[2][kintold] ) |
| 543 |
#print nodecheck, Node[0][nodecheck],Node[1][nodecheck],Node[2][nodecheck] |
| 544 |
#nodecheck=int(Element_Nodes[3][kintfault] ) |
| 545 |
#print nodecheck, Node[0][nodecheck],Node[1][nodecheck],Node[2][nodecheck] |
| 546 |
|
| 547 |
#nodecheck=int(Element_Nodes[6][kintold] ) |
| 548 |
#print nodecheck, Node[0][nodecheck],Node[1][nodecheck],Node[2][nodecheck] |
| 549 |
#nodecheck=int(Element_Nodes[7][kintfault] ) |
| 550 |
#print nodecheck, Node[0][nodecheck],Node[1][nodecheck],Node[2][nodecheck] |
| 551 |
|
| 552 |
#nodecheck=int(Element_Nodes[5][kintold] ) |
| 553 |
#print nodecheck, Node[0][nodecheck],Node[1][nodecheck],Node[2][nodecheck] |
| 554 |
#nodecheck=int(Element_Nodes[4][kintfault] ) |
| 555 |
#print nodecheck, Node[0][nodecheck],Node[1][nodecheck],Node[2][nodecheck] |
| 556 |
|
| 557 |
ContactElement_ref[k]=k |
| 558 |
ContactElement_tag[k]=2 |
| 559 |
|
| 560 |
ContactElement_Nodes[0][k]=Element_Nodes[1][kintold] |
| 561 |
ContactElement_Nodes[1][k]=Element_Nodes[2][kintold] |
| 562 |
ContactElement_Nodes[2][k]=Element_Nodes[6][kintold] |
| 563 |
ContactElement_Nodes[3][k]=Element_Nodes[5][kintold] |
| 564 |
ContactElement_Nodes[4][k]=Element_Nodes[0][kintold] |
| 565 |
ContactElement_Nodes[5][k]=Element_Nodes[3][kintold] |
| 566 |
ContactElement_Nodes[6][k]=Element_Nodes[7][kintold] |
| 567 |
ContactElement_Nodes[7][k]=Element_Nodes[4][kintold] |
| 568 |
|
| 569 |
ContactElement_Nodes[8][k]=Element_Nodes[0][kintfault] |
| 570 |
ContactElement_Nodes[9][k]=Element_Nodes[3][kintfault] |
| 571 |
ContactElement_Nodes[10][k]=Element_Nodes[7][kintfault] |
| 572 |
ContactElement_Nodes[11][k]=Element_Nodes[4][kintfault] |
| 573 |
ContactElement_Nodes[12][k]=Element_Nodes[1][kintfault] |
| 574 |
ContactElement_Nodes[13][k]=Element_Nodes[2][kintfault] |
| 575 |
ContactElement_Nodes[14][k]=Element_Nodes[6][kintfault] |
| 576 |
ContactElement_Nodes[15][k]=Element_Nodes[5][kintfault] |
| 577 |
|
| 578 |
meshfaultL.write("%s %d\n"%(ContactElement_Type,ContactElement_Num)) |
| 579 |
|
| 580 |
for i in range(ContactElement_Num): |
| 581 |
meshfaultL.write("%d %d"%(ContactElement_ref[i],ContactElement_tag[i])) |
| 582 |
for j in range(ContactElement_numNodes): |
| 583 |
meshfaultL.write(" %d"%ContactElement_Nodes[j][i]) |
| 584 |
meshfaultL.write("\n") |
| 585 |
|
| 586 |
# point sources (not supported yet) |
| 587 |
|
| 588 |
meshfaultL.write("Point1 0") |
| 589 |
|
| 590 |
|
| 591 |
|
| 592 |
meshfaultL.close() |
| 593 |
|
| 594 |
mydomainfile = faultL(width,width, height,ne, ne, ne,contact=True,xstart=fstart,xend=fend) |
| 595 |
|
Parent Directory
| 
Revision Log