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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 884 - (hide annotations)
Mon Oct 30 06:37:30 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: 22234 byte(s)
simple mesh generator for slip_stress
1 gross 884 # $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