/[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 887 - (show annotations)
Thu Nov 2 07:17:07 2006 UTC (13 years, 4 months ago) by gross
Original Path: trunk/finley/test/python/slip_stress_mesh.py
File MIME type: text/x-python
File size: 22412 byte(s)
some more work on the slip distribution
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

  ViewVC Help
Powered by ViewVC 1.1.26