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

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

Parent Directory Parent Directory | Revision Log Revision Log


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

  ViewVC Help
Powered by ViewVC 1.1.26