/[escript]/trunk/finley/test/python/slip_stress_mesh_old.py
ViewVC logotype

Contents of /trunk/finley/test/python/slip_stress_mesh_old.py

Parent Directory Parent Directory | Revision Log Revision Log


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

  ViewVC Help
Powered by ViewVC 1.1.26