/[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 888 - (hide annotations)
Tue Nov 7 08:31:26 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: 22453 byte(s)
Problem in Mesh_findMatchingFaces.c fixed:
default tolerance in python interface was set too tight.



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 gross 887 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 gross 884
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 gross 887 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 gross 884 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 gross 888 # ne_w=int((ne/height)*width+0.5)
595     ne_w=ne
596 gross 887 mydomainfile = faultL(width,width, height,ne, ne, ne,contact=True,xstart=fstart,xend=fend)

  ViewVC Help
Powered by ViewVC 1.1.26