/[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 893 - (hide annotations)
Wed Nov 8 08:20:19 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: 24698 byte(s)
small bug fixed
1 gross 884 # $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 gross 893 ne = 10
33 gross 887 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 gross 884
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 gross 893
119 gross 884 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 gross 893 k= i0 + N0*i1 + N0*N1*i2 # M0*i0+M1*i1+M2*i2;
141     Node_ref[k]= i0 + N0*i1 + N0*N1*i2
142 gross 884 # 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 gross 893 Node[2][k]=(i2)*l2/(N2-1)
148 gross 884
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 gross 893 # CHANGED:
165     k=Fault_NE+i1+(n1double-1)*i2
166     Node_ref[k]= k #Fault_NE + i1 + (n1double-1)*i2
167 gross 884 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 gross 893 Node[2][k]=i2start*l2/ne2 + (i2+1)*l2/ne2
172 gross 884 # elif n1double==0:
173     # elif n2double==0:
174 gross 893 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 gross 884
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 gross 893 print "M0,M1,M2",M0,M1,M2
197    
198 gross 884 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 gross 893 #CHANGED:
208 gross 884 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 gross 893 #print "x0s,x1s,x2s,x0e,x1e,x2e", x0s,x1s,x2s,x0e,x1e,x2e
226 gross 884 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 gross 893 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 gross 884 else:
236     kfault=0.
237     # determine bottom corner node of each element
238    
239     # Ref for normal interior node:
240 gross 893 kold=int(i0start+N0*(i1start + i1) + (N0*N1)*(i2start+i2))
241 gross 884 #print kold, Node[0][kold],Node[1][kold],Node[2][kold]
242     # Ref for interior element:
243 gross 893 kint=int(i0start + ne0*(i1start+i1) + (ne0*ne1)*(i2start+i2))
244 gross 884 #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 gross 893 kfaultref=int(Fault_NE+i1+(n1double-1)*i2)
258 gross 884 elif (i1==0):
259     # nearest fault node
260 gross 893 kfaultref=int(Fault_NE+i1+(i2-1.)*(n1double-1))
261 gross 884 elif (i2==0):
262     # nearest fault node
263 gross 893 kfaultref=int(Fault_NE+(i1-1.) + i2*(n1double-1))
264 gross 884 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 gross 893 #nodecheck=int(Element_Nodes[7][kint] )
275     #print nodecheck, Node[0][nodecheck],Node[1][nodecheck],Node[2][nodecheck]
276    
277 gross 884 Element_Nodes[7][kint]=kfaultref
278    
279 gross 893 #print kfaultref, Node[0][kfaultref],Node[1][kfaultref],Node[2][kfaultref]
280    
281    
282 gross 884 elif (i1==0 and i2==n2double-1):
283 gross 893
284     #nodecheck=int(Element_Nodes[3][kint] )
285     #print nodecheck, Node[0][nodecheck],Node[1][nodecheck],Node[2][nodecheck]
286    
287 gross 884 Element_Nodes[3][kint]=kfaultref
288    
289 gross 893 #print kfaultref, Node[0][kfaultref],Node[1][kfaultref],Node[2][kfaultref]
290    
291 gross 884 elif (i1==n1double-1 and i2==0):
292 gross 893 #nodecheck=int(Element_Nodes[4][kint] )
293     #print nodecheck, Node[0][nodecheck],Node[1][nodecheck],Node[2][nodecheck]
294    
295 gross 884 Element_Nodes[4][kint]=kfaultref
296 gross 893 #print kfaultref, Node[0][kfaultref],Node[1][kfaultref],Node[2][kfaultref]
297 gross 884
298     elif (i1==n1double-1 and i2==n2double-1):
299 gross 893 #nodecheck=int(Element_Nodes[0][kint] )
300     #print nodecheck, Node[0][nodecheck],Node[1][nodecheck],Node[2][nodecheck]
301    
302 gross 884 Element_Nodes[0][kint]=kfaultref
303 gross 893 #print kfaultref, Node[0][kfaultref],Node[1][kfaultref],Node[2][kfaultref]
304 gross 884
305     # overwrite 4 sides of fault (only 2 nodes changed)
306     elif (i1==0):
307 gross 893
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 gross 884 Element_Nodes[3][kint]=kfaultref
314 gross 893 kfaultref1=int(kfaultref+(n1double-1))
315 gross 884 Element_Nodes[7][kint]=kfaultref1
316    
317 gross 893 #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 gross 884 elif (i1==n1double-1):
323 gross 893
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 gross 884 Element_Nodes[0][kint]=kfaultref
331 gross 893 kfaultref1=kfaultref+(n1double-1)
332 gross 884 Element_Nodes[4][kint]=kfaultref1
333    
334 gross 893 #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 gross 884 elif (i2==0):
338 gross 893
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 gross 884 Element_Nodes[4][kint]=kfaultref
345 gross 893 kfaultref1=kfaultref+1
346 gross 884 Element_Nodes[7][kint]=kfaultref1
347    
348 gross 893 #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 gross 884 elif (i2==n2double-1):
352 gross 893
353 gross 884 #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 gross 893 Element_Nodes[0][kint]=kfaultref
359     kfaultref1=kfaultref+1
360 gross 884 Element_Nodes[3][kint]=kfaultref1
361    
362 gross 893 #print kfaultref, Node[0][kfaultref],Node[1][kfaultref],Node[2][kfaultref]
363     #print kfaultref1, Node[0][kfaultref1],Node[1][kfaultref1],Node[2][kfaultref1]
364 gross 884
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 gross 893 kfaultref1=kfaultref+1
374     kfaultref2=kfaultref+(n1double-1)
375     kfaultref3=kfaultref+1+(n1double-1)
376 gross 884
377 gross 893 #nodecheck=int(Element_Nodes[3][kint] )
378 gross 884 #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 gross 893 #nodecheck=int(Element_Nodes[4][kint] )
382 gross 884 #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 gross 893 #nodecheck=int(Element_Nodes[7][kint] )
386 gross 884 #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 gross 893 kintold=int((i0start-1) + ne0*(i1start+i1) + ne0*ne1*(i2start+i2))
576 gross 884
577     # define reference for interior elements with x0>x0s
578     # here the double nodes are the fault nodes
579    
580 gross 893 kintfault=int(i0start + ne0*(i1start+i1) + ne0*ne1*(i2start+i2))
581 gross 884
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 gross 893
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