/[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 1312 - (hide annotations)
Mon Sep 24 06:18:44 2007 UTC (12 years, 5 months ago) by ksteube
File MIME type: text/x-python
File size: 25119 byte(s)
The MPI branch is hereby closed. All future work should be in trunk.

Previously in revision 1295 I merged the latest changes to trunk into trunk-mpi-branch.
In this revision I copied all files from trunk-mpi-branch over the corresponding
trunk files. I did not use 'svn merge', it was a copy.

1 ksteube 1312 #
2 gross 884 # $Id$
3 ksteube 1312 #
4     #######################################################
5     #
6     # Copyright 2003-2007 by ACceSS MNRF
7     # Copyright 2007 by University of Queensland
8     #
9     # http://esscc.uq.edu.au
10     # Primary Business: Queensland, Australia
11     # Licensed under the Open Software License version 3.0
12     # http://www.opensource.org/licenses/osl-3.0.php
13     #
14     #######################################################
15     #
16    
17 gross 884 """
18    
19     generates finley mesh simple vertical fault
20    
21     THIS CODE CREATES RICH CONTACT ELEMENTS AND RICH FACE ELEMENTS
22     with fix for contact elements at FAULT ENDS
23    
24    
25     @var __author__: name of author
26     @var __copyright__: copyrights
27     @var __license__: licence agreement
28     @var __url__: url entry point on documentation
29     @var __version__: version
30     @var __date__: date of the version
31     """
32    
33     __author__="Louise Kettle"
34     __copyright__=""" Copyright (c) 2006 by ACcESS MNRF
35     http://www.access.edu.au
36     Primary Business: Queensland, Australia"""
37     __license__="""Licensed under the Open Software License version 3.0
38     http://www.opensource.org/licenses/osl-3.0.php"""
39     __url__="http://www.iservo.edu.au/esys"
40     __version__="$Revision$"
41     __date__="$Date$"
42    
43     from esys.escript import *
44     from numarray import zeros,Float,array,size
45    
46     #... generate domain ...
47 gross 893 ne = 10
48 gross 887 width = 100000.
49     height = 30000.
50     fstart=numarray.array([width/2.,7.*width/16.,3.*height/8.])
51     fend=numarray.array([width/2.,9.*width/16.,5.*height/8.])
52 gross 884
53     def faultL(l0,l1, l2,ne0, ne1, ne2,contact=False,xstart=zeros(3),xend=zeros(3)):
54     meshfaultL=open('meshfault3D.fly','w')
55    
56     FaultError1="ERROR: fault defined on or too close to an outer surface"
57     FaultError2="ERROR: the mesh is too coarse for fault"
58    
59     N0=ne0+1
60     N1=ne1+1
61     N2=ne2+1
62     if (N0<=N1 and N0<=N2):
63     if (N1 <= N2):
64     M0=1
65     M1=N0
66     M2=N0*N1
67     M0i=1
68     M1i=ne0
69     M2i=ne0*ne1
70     else:
71     M0=1
72     M2=N0
73     M1=N0*N2
74     M0i=1
75     M2i=ne0
76     M1i=ne0*ne2
77     elif (N1<=N2 and N1<=N0):
78     if (N2 <= N0):
79     M1=1
80     M2=N1
81     M0=N2*N1
82     M1i=1
83     M2i=ne1
84     M0i=ne2*ne1
85     else:
86     M1=1
87     M0=N1
88     M2=N1*N0
89     M1i=1
90     M0i=ne1
91     M2i=ne0*ne1
92     else:
93     if (N0 <= N1):
94     M2=1
95     M0=N2
96     M1=N2*N0
97     M2i=1
98     M0i=ne2
99     M1i=ne0*ne2
100     else:
101     M2=1
102     M1=N2
103     M0=N1*N2
104     M2i=1
105     M1i=ne2
106     M0i=ne2*ne1
107    
108     dim=3
109     Element_numNodes=8
110     Element_Num=ne0*ne1*ne2
111     if contact==False:
112     numNodes=N0*N1*N2
113    
114     else:
115     # define double (contact element) nodes on interior of fault
116     i0start=round(xstart[0]*ne0/l0)
117     i1start=round(xstart[1]*ne1/l1)
118     i2start=round(xstart[2]*ne2/l2)
119     i0end=round(xend[0]*ne0/l0)
120     i1end=round(xend[1]*ne1/l1)
121     i2end=round(xend[2]*ne2/l2)
122     n0double=int(i0end)-int(i0start)
123     n1double=int(i1end)-int(i1start)
124     n2double=int(i2end)-int(i2start)
125     if (i0start == 0) or (i1start==0) or (i2start==0):
126     raise FaultError1
127    
128     if (i0end == ne0) or (i1end==ne1) or (i2end==ne2):
129     raise FaultError1
130    
131     if n0double==0:
132     numNodes=N0*N1*N2+(n1double-1)*(n2double-1)
133 gross 893
134 gross 884 elif n1double==0:
135     numNodes=N0*N1*N2+(n0double-1)*(n2double-1)
136    
137     elif n2double==0:
138     numNodes=N0*N1*N2+(n0double-1)*(n1double-1)
139    
140     # define nodes for normal elements
141     # there are N0*N1*N2 normal nodes
142    
143     Node=zeros([3,numNodes],Float)
144     Node_ref=zeros(numNodes,Float)
145     Node_DOF=zeros(numNodes,Float)
146     Node_tag=zeros(numNodes,Float)
147    
148     meshfaultL.write("KettleFault\n")
149     #print 'Nodes'
150     meshfaultL.write("%dD-nodes %d\n"%(dim,numNodes))
151    
152     for i2 in range(N2):
153     for i1 in range (N1):
154     for i0 in range(N0):
155 gross 893 k= i0 + N0*i1 + N0*N1*i2 # M0*i0+M1*i1+M2*i2;
156     Node_ref[k]= i0 + N0*i1 + N0*N1*i2
157 gross 884 # no periodic boundary conditions
158     Node_DOF[k]=Node_ref[k]
159     Node_tag[k]=0
160     Node[0][k]=(i0)*l0/(N0-1)
161     Node[1][k]=(i1)*l1/(N1-1)
162 gross 893 Node[2][k]=(i2)*l2/(N2-1)
163 gross 884
164     # define double nodes on fault (will have same coordinates as some of nodes already defined)
165     # only get double nodes on INTERIOR of fault
166    
167     if contact==True:
168     Fault_NE=N0*N1*N2
169     if n0double==0:
170     if(n1double<=n2double):
171     M1f=1
172     M2f=n1double-1
173     else:
174     M2f=1
175     M1f=n2double-1
176    
177     for i2 in range(n2double-1):
178     for i1 in range(n1double-1):
179 gross 893 # CHANGED:
180     k=Fault_NE+i1+(n1double-1)*i2
181     Node_ref[k]= k #Fault_NE + i1 + (n1double-1)*i2
182 gross 884 Node_DOF[k]=Node_ref[k]
183     Node_tag[k]=1
184     Node[0][k]=i0start*l0/ne0
185     Node[1][k]=i1start*l1/ne1 + (i1+1)*l1/ne1
186 gross 893 Node[2][k]=i2start*l2/ne2 + (i2+1)*l2/ne2
187 gross 884 # elif n1double==0:
188     # elif n2double==0:
189 gross 893 print "fstart = ",[i0start*l0/ne0, i1start*l1/ne1 , i2start*l2/ne2]
190     print "fend = ", [i0start*l0/ne0 , i1start*l1/ne1 + n1double*l1/ne1, i2start*l2/ne2 + n2double*l2/ne2]
191 gross 884
192     # write nodes to file
193     for i in range(numNodes):
194     meshfaultL.write("%d %d %d"%(Node_ref[i],Node_DOF[i],Node_tag[i]))
195     for j in range(dim):
196     meshfaultL.write(" %lf"%Node[j][i])
197     meshfaultL.write("\n")
198    
199    
200    
201    
202     # defining interior elements
203     # there are ne0*ne1*ne2 interior elements
204    
205     Element_Nodes=zeros([8,ne0*ne1*ne2],Float)
206     Element_ref=zeros(ne0*ne1*ne2,Float)
207     Element_tag=zeros(ne0*ne1*ne2,Float)
208    
209     #print 'Interior elements'
210    
211 gross 893 print "M0,M1,M2",M0,M1,M2
212    
213 gross 884 for i2 in range(ne2):
214     for i1 in range (ne1):
215     for i0 in range(ne0):
216     k=i0 + ne0*i1 + ne0*ne1*i2;
217     # define corner node (node0)
218     node0=i0 + N0*i1 + N0*N1*i2;
219     Element_ref[k]=k
220     Element_tag[k]=0
221     # for hex8 the interior elements are specified by 8 nodes
222 gross 893 #CHANGED:
223 gross 884 Element_Nodes[0][k]=node0;
224     Element_Nodes[1][k]=node0+1;
225     Element_Nodes[2][k]=node0+N0+1;
226     Element_Nodes[3][k]=node0+N0;
227     Element_Nodes[4][k]=node0+N0*N1;
228     Element_Nodes[5][k]=node0+N0*N1+1;
229     Element_Nodes[6][k]=node0+N0*N1+N0+1;
230     Element_Nodes[7][k]=node0+N0*N1+N0;
231    
232     if contact==True:
233     if n0double==0:
234     x0s= i0start*l0/ne0
235     x1s= i1start*l1/ne1
236     x2s= i2start*l2/ne2
237     x0e= i0end*l0/ne0
238     x1e= i1end*l1/ne1
239     x2e= i2end*l2/ne2
240 gross 893 #print "x0s,x1s,x2s,x0e,x1e,x2e", x0s,x1s,x2s,x0e,x1e,x2e
241 gross 884 if (n1double==1) or (n2double==1):
242     raise FaultError2
243     for i2 in range(n2double):
244     for i1 in range(n1double):
245     # here the coordinates of kfault and kold are the same
246     # Ref for fault node (only on interior nodes of fault):
247     if (i1>0) and (i2>0):
248 gross 893 kfault=Fault_NE+(i1-1.)+(n1double-1)*(i2-1.)
249     #print kfault , Node[0][int(kfault)],Node[1][int(kfault)],Node[2][int(kfault)]
250 gross 884 else:
251     kfault=0.
252     # determine bottom corner node of each element
253    
254     # Ref for normal interior node:
255 gross 893 kold=int(i0start+N0*(i1start + i1) + (N0*N1)*(i2start+i2))
256 gross 884 #print kold, Node[0][kold],Node[1][kold],Node[2][kold]
257     # Ref for interior element:
258 gross 893 kint=int(i0start + ne0*(i1start+i1) + (ne0*ne1)*(i2start+i2))
259 gross 884 #print kint, Element_Nodes[0][kint]
260    
261     x0= (i0start)*l0/ne0
262     x1= (i1start+i1)*l1/ne1
263     x2= (i2start+i2)*l2/ne2
264    
265     # for x0 > xstart we need to overwrite old Nodes in interior element references
266     # with fault nodes:
267    
268     # for the interior elements with x1<x1s and x2<x2s the only nodes need changing
269     # are on the fault:
270     if (i1==0) and (i2==0):
271     # nearest fault node:
272 gross 893 kfaultref=int(Fault_NE+i1+(n1double-1)*i2)
273 gross 884 elif (i1==0):
274     # nearest fault node
275 gross 893 kfaultref=int(Fault_NE+i1+(i2-1.)*(n1double-1))
276 gross 884 elif (i2==0):
277     # nearest fault node
278 gross 893 kfaultref=int(Fault_NE+(i1-1.) + i2*(n1double-1))
279 gross 884 else:
280     # looking at element with fault node on bottom corner
281     kfaultref=int(kfault)
282    
283     #print x0,x1,x2
284     #print kold, Node[0][kold],Node[1][kold],Node[2][kold]
285     #print kfaultref, Node[0][kfaultref],Node[1][kfaultref],Node[2][kfaultref]
286    
287     # overwrite 4 outer corner elements of fault (only one node changed)
288     if (i1==0 and i2==0):
289 gross 893 #nodecheck=int(Element_Nodes[7][kint] )
290     #print nodecheck, Node[0][nodecheck],Node[1][nodecheck],Node[2][nodecheck]
291    
292 gross 884 Element_Nodes[7][kint]=kfaultref
293    
294 gross 893 #print kfaultref, Node[0][kfaultref],Node[1][kfaultref],Node[2][kfaultref]
295    
296    
297 gross 884 elif (i1==0 and i2==n2double-1):
298 gross 893
299     #nodecheck=int(Element_Nodes[3][kint] )
300     #print nodecheck, Node[0][nodecheck],Node[1][nodecheck],Node[2][nodecheck]
301    
302 gross 884 Element_Nodes[3][kint]=kfaultref
303    
304 gross 893 #print kfaultref, Node[0][kfaultref],Node[1][kfaultref],Node[2][kfaultref]
305    
306 gross 884 elif (i1==n1double-1 and i2==0):
307 gross 893 #nodecheck=int(Element_Nodes[4][kint] )
308     #print nodecheck, Node[0][nodecheck],Node[1][nodecheck],Node[2][nodecheck]
309    
310 gross 884 Element_Nodes[4][kint]=kfaultref
311 gross 893 #print kfaultref, Node[0][kfaultref],Node[1][kfaultref],Node[2][kfaultref]
312 gross 884
313     elif (i1==n1double-1 and i2==n2double-1):
314 gross 893 #nodecheck=int(Element_Nodes[0][kint] )
315     #print nodecheck, Node[0][nodecheck],Node[1][nodecheck],Node[2][nodecheck]
316    
317 gross 884 Element_Nodes[0][kint]=kfaultref
318 gross 893 #print kfaultref, Node[0][kfaultref],Node[1][kfaultref],Node[2][kfaultref]
319 gross 884
320     # overwrite 4 sides of fault (only 2 nodes changed)
321     elif (i1==0):
322 gross 893
323     #nodecheck=int(Element_Nodes[3][kint] )
324     #print nodecheck, Node[0][nodecheck],Node[1][nodecheck],Node[2][nodecheck]
325     #nodecheck=int(Element_Nodes[7][kint] )
326     #print nodecheck, Node[0][nodecheck],Node[1][nodecheck],Node[2][nodecheck]
327    
328 gross 884 Element_Nodes[3][kint]=kfaultref
329 gross 893 kfaultref1=int(kfaultref+(n1double-1))
330 gross 884 Element_Nodes[7][kint]=kfaultref1
331    
332 gross 893 #print kfaultref, Node[0][kfaultref],Node[1][kfaultref],Node[2][kfaultref]
333     #print kfaultref1, Node[0][kfaultref1],Node[1][kfaultref1],Node[2][kfaultref1]
334    
335    
336    
337 gross 884 elif (i1==n1double-1):
338 gross 893
339     #nodecheck=int(Element_Nodes[0][kint] )
340     #print nodecheck, Node[0][nodecheck],Node[1][nodecheck],Node[2][nodecheck]
341     #nodecheck=int(Element_Nodes[4][kint] )
342     #print nodecheck, Node[0][nodecheck],Node[1][nodecheck],Node[2][nodecheck]
343    
344    
345 gross 884 Element_Nodes[0][kint]=kfaultref
346 gross 893 kfaultref1=kfaultref+(n1double-1)
347 gross 884 Element_Nodes[4][kint]=kfaultref1
348    
349 gross 893 #print kfaultref, Node[0][kfaultref],Node[1][kfaultref],Node[2][kfaultref]
350     #print kfaultref1, Node[0][kfaultref1],Node[1][kfaultref1],Node[2][kfaultref1]
351    
352 gross 884 elif (i2==0):
353 gross 893
354     #nodecheck=int(Element_Nodes[4][kint] )
355     #print nodecheck, Node[0][nodecheck],Node[1][nodecheck],Node[2][nodecheck]
356     #nodecheck=int(Element_Nodes[7][kint] )
357     #print nodecheck, Node[0][nodecheck],Node[1][nodecheck],Node[2][nodecheck]
358    
359 gross 884 Element_Nodes[4][kint]=kfaultref
360 gross 893 kfaultref1=kfaultref+1
361 gross 884 Element_Nodes[7][kint]=kfaultref1
362    
363 gross 893 #print kfaultref, Node[0][kfaultref],Node[1][kfaultref],Node[2][kfaultref]
364     #print kfaultref1, Node[0][kfaultref1],Node[1][kfaultref1],Node[2][kfaultref1]
365    
366 gross 884 elif (i2==n2double-1):
367 gross 893
368 gross 884 #nodecheck=int(Element_Nodes[0][kint] )
369     #print nodecheck, Node[0][nodecheck],Node[1][nodecheck],Node[2][nodecheck]
370     #nodecheck=int(Element_Nodes[3][kint] )
371     #print nodecheck, Node[0][nodecheck],Node[1][nodecheck],Node[2][nodecheck]
372    
373 gross 893 Element_Nodes[0][kint]=kfaultref
374     kfaultref1=kfaultref+1
375 gross 884 Element_Nodes[3][kint]=kfaultref1
376    
377 gross 893 #print kfaultref, Node[0][kfaultref],Node[1][kfaultref],Node[2][kfaultref]
378     #print kfaultref1, Node[0][kfaultref1],Node[1][kfaultref1],Node[2][kfaultref1]
379 gross 884
380     # overwrite interior fault elements (4 nodes changed)
381     else:
382     #nodecheck=int(Element_Nodes[0][kint] )
383     #print nodecheck, Node[0][nodecheck],Node[1][nodecheck],Node[2][nodecheck]
384     #print kfaultref, Node[0][kfaultref],Node[1][kfaultref],Node[2][kfaultref]
385    
386     Element_Nodes[0][kint]=kfaultref
387     #if (x1<x1e and x2<x2e):
388 gross 893 kfaultref1=kfaultref+1
389     kfaultref2=kfaultref+(n1double-1)
390     kfaultref3=kfaultref+1+(n1double-1)
391 gross 884
392 gross 893 #nodecheck=int(Element_Nodes[3][kint] )
393 gross 884 #print nodecheck, Node[0][nodecheck],Node[1][nodecheck],Node[2][nodecheck]
394     #print kfaultref1, Node[0][kfaultref1],Node[1][kfaultref1],Node[2][kfaultref1]
395    
396 gross 893 #nodecheck=int(Element_Nodes[4][kint] )
397 gross 884 #print nodecheck, Node[0][nodecheck],Node[1][nodecheck],Node[2][nodecheck]
398     #print kfaultref2, Node[0][kfaultref2],Node[1][kfaultref2],Node[2][kfaultref2]
399    
400 gross 893 #nodecheck=int(Element_Nodes[7][kint] )
401 gross 884 #print nodecheck, Node[0][nodecheck],Node[1][nodecheck],Node[2][nodecheck]
402     #print kfaultref3, Node[0][kfaultref3],Node[1][kfaultref3],Node[2][kfaultref3]
403    
404    
405     Element_Nodes[3][kint]=kfaultref1
406     Element_Nodes[4][kint]=kfaultref2
407     Element_Nodes[7][kint]=kfaultref3
408     #elif x1<x1e:
409     # kfaultref1=kfaultref+M1f
410     # Element_Nodes[3][kint]=kfaultref1
411     #elif x2<x2e:
412     # kfaultref2=kfaultref+M2f
413     # Element_Nodes[4][kint]=kfaultref2
414    
415    
416     # print kint, kfaultref
417    
418     # write interior elements to file
419     Element_Type='Hex8'
420     meshfaultL.write("%s %d\n"%(Element_Type,Element_Num))
421    
422     for i in range(Element_Num):
423     meshfaultL.write("%d %d"%(Element_ref[i],Element_tag[i]))
424     for j in range(Element_numNodes):
425     meshfaultL.write(" %d"%Element_Nodes[j][i])
426     meshfaultL.write("\n")
427    
428     # face elements
429     FaceElement_Type='Hex8Face'
430     FaceElement_Num= 2*(ne0*ne1 + ne0*ne2 + ne1*ne2)
431     FaceElement_numNodes=8
432    
433     meshfaultL.write("%s %d\n"%(FaceElement_Type,FaceElement_Num))
434    
435     FaceElement_Nodes=zeros([FaceElement_numNodes,FaceElement_Num],Float)
436     FaceElement_ref=zeros(FaceElement_Num,Float)
437     FaceElement_tag=zeros(FaceElement_Num,Float)
438    
439     kcount=0
440    
441     # defining face elements on x2=0 face
442     for i1 in range (ne1):
443     for i0 in range(ne0):
444     i2=0
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;
452     FaceElement_Nodes[1][kcount]=node0+1;
453     FaceElement_Nodes[2][kcount]=node0+N0+1;
454     FaceElement_Nodes[3][kcount]=node0+N0;
455     FaceElement_Nodes[4][kcount]=node0+N0*N1;
456     FaceElement_Nodes[5][kcount]=node0+N0*N1+1;
457     FaceElement_Nodes[6][kcount]=node0+N0*N1+N0+1;
458     FaceElement_Nodes[7][kcount]=node0+N0*N1+N0;
459     kcount+=1
460    
461     # defining face elements on x2=L face
462     for i1 in range (ne1):
463     for i0 in range(ne0):
464     i2=ne2-1
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+N0*N1;
472     FaceElement_Nodes[1][kcount]=node0+N0*N1+1;
473     FaceElement_Nodes[2][kcount]=node0+N0*N1+N0+1;
474     FaceElement_Nodes[3][kcount]=node0+N0*N1+N0;
475     FaceElement_Nodes[4][kcount]=node0;
476     FaceElement_Nodes[5][kcount]=node0+1;
477     FaceElement_Nodes[6][kcount]=node0+N0+1;
478     FaceElement_Nodes[7][kcount]=node0+N0;
479     kcount+=1
480    
481     # defining face elements on x1=0 face
482     for i2 in range (ne2):
483     for i0 in range(ne0):
484     i1=0
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;
492     FaceElement_Nodes[1][kcount]=node0+N0*N1;
493     FaceElement_Nodes[2][kcount]=node0+N0*N1+1;
494     FaceElement_Nodes[3][kcount]=node0+1;
495     FaceElement_Nodes[4][kcount]=node0+N0;
496     FaceElement_Nodes[5][kcount]=node0+N0*N1+N0;
497     FaceElement_Nodes[6][kcount]=node0+N0*N1+N0+1;
498     FaceElement_Nodes[7][kcount]=node0+N0+1;
499     kcount+=1
500    
501     # defining face elements on x1=L face
502     for i2 in range (ne2):
503     for i0 in range(ne0):
504     i1=ne1-1
505     k=i0 + ne0*i1 + ne0*ne1*i2;
506     # define corner node (node0)
507     node0=i0 + N0*i1 + N0*N1*i2;
508     FaceElement_ref[kcount]=kcount
509     FaceElement_tag[kcount]=3
510     # for hex8face the face elements are specified by 8 nodes
511     FaceElement_Nodes[0][kcount]=node0+N0;
512     FaceElement_Nodes[1][kcount]=node0+N0*N1+N0;
513     FaceElement_Nodes[2][kcount]=node0+N0*N1+N0+1;
514     FaceElement_Nodes[3][kcount]=node0+N0+1;
515     FaceElement_Nodes[4][kcount]=node0;
516     FaceElement_Nodes[5][kcount]=node0+N0*N1;
517     FaceElement_Nodes[6][kcount]=node0+N0*N1+1;
518     FaceElement_Nodes[7][kcount]=node0+1;
519     kcount+=1
520    
521     # defining face elements on x0=0 face
522     for i2 in range (ne2):
523     for i1 in range(ne1):
524     i0=0
525     k=i0 + ne0*i1 + ne0*ne1*i2;
526     # define corner node (node0)
527     node0=i0 + N0*i1 + N0*N1*i2;
528     FaceElement_ref[kcount]=kcount
529     FaceElement_tag[kcount]=3
530     # for hex8face the face elements are specified by 8 nodes
531     FaceElement_Nodes[0][kcount]=node0;
532     FaceElement_Nodes[1][kcount]=node0+N0;
533     FaceElement_Nodes[2][kcount]=node0+N0*N1+N0;
534     FaceElement_Nodes[3][kcount]=node0+N0*N1;
535     FaceElement_Nodes[4][kcount]=node0+1;
536     FaceElement_Nodes[5][kcount]=node0+N0+1;
537     FaceElement_Nodes[6][kcount]=node0+N0*N1+N0+1;
538     FaceElement_Nodes[7][kcount]=node0+N0*N1+1;
539     kcount+=1
540    
541     # defining face elements on x0=L face
542     for i2 in range (ne2):
543     for i1 in range(ne1):
544     i0=ne1-1
545     k=i0 + ne0*i1 + ne0*ne1*i2;
546     # define corner node (node0)
547     node0=i0 + N0*i1 + N0*N1*i2;
548     FaceElement_ref[kcount]=kcount
549     FaceElement_tag[kcount]=3
550     # for hex8face the face elements are specified by 8 nodes
551     FaceElement_Nodes[0][kcount]=node0+1;
552     FaceElement_Nodes[1][kcount]=node0+N0+1;
553     FaceElement_Nodes[2][kcount]=node0+N0*N1+N0+1;
554     FaceElement_Nodes[3][kcount]=node0+N0*N1+1;
555     FaceElement_Nodes[4][kcount]=node0;
556     FaceElement_Nodes[5][kcount]=node0+N0;
557     FaceElement_Nodes[6][kcount]=node0+N0*N1+N0;
558     FaceElement_Nodes[7][kcount]=node0+N0*N1;
559     kcount+=1
560    
561    
562    
563    
564     for i in range(FaceElement_Num):
565     meshfaultL.write("%d %d"%(FaceElement_ref[i],FaceElement_tag[i]))
566     for j in range(FaceElement_numNodes):
567     meshfaultL.write(" %d"%FaceElement_Nodes[j][i])
568     meshfaultL.write("\n")
569    
570    
571     # contact elements
572     ContactElement_Type='Hex8Face_Contact'
573     ContactElement_Num=0
574     ContactElement_numNodes=16
575     # print contact elements on fault
576     if contact==True:
577     if n0double==0:
578     ContactElement_Num=(n1double)*(n2double)
579     ContactElement_Nodes=zeros([ContactElement_numNodes,ContactElement_Num],Float)
580     ContactElement_ref=zeros(ContactElement_Num,Float)
581     ContactElement_tag=zeros(ContactElement_Num,Float)
582     #print ContactElement_Num
583    
584     for i2 in range(n2double):
585     for i1 in range(n1double):
586     k=i1+(n1double)*i2
587     #print k
588     # define reference for interior elements with x0<=x0s
589     # here the nodes are the old interior nodes
590 gross 893 kintold=int((i0start-1) + ne0*(i1start+i1) + ne0*ne1*(i2start+i2))
591 gross 884
592     # define reference for interior elements with x0>x0s
593     # here the double nodes are the fault nodes
594    
595 gross 893 kintfault=int(i0start + ne0*(i1start+i1) + ne0*ne1*(i2start+i2))
596 gross 884
597     #nodecheck=int(Element_Nodes[1][kintold] )
598     #print nodecheck, Node[0][nodecheck],Node[1][nodecheck],Node[2][nodecheck]
599     #nodecheck=int(Element_Nodes[0][kintfault] )
600     #print nodecheck, Node[0][nodecheck],Node[1][nodecheck],Node[2][nodecheck]
601    
602     #nodecheck=int(Element_Nodes[2][kintold] )
603     #print nodecheck, Node[0][nodecheck],Node[1][nodecheck],Node[2][nodecheck]
604     #nodecheck=int(Element_Nodes[3][kintfault] )
605     #print nodecheck, Node[0][nodecheck],Node[1][nodecheck],Node[2][nodecheck]
606    
607     #nodecheck=int(Element_Nodes[6][kintold] )
608     #print nodecheck, Node[0][nodecheck],Node[1][nodecheck],Node[2][nodecheck]
609     #nodecheck=int(Element_Nodes[7][kintfault] )
610     #print nodecheck, Node[0][nodecheck],Node[1][nodecheck],Node[2][nodecheck]
611    
612     #nodecheck=int(Element_Nodes[5][kintold] )
613     #print nodecheck, Node[0][nodecheck],Node[1][nodecheck],Node[2][nodecheck]
614     #nodecheck=int(Element_Nodes[4][kintfault] )
615     #print nodecheck, Node[0][nodecheck],Node[1][nodecheck],Node[2][nodecheck]
616    
617     ContactElement_ref[k]=k
618     ContactElement_tag[k]=2
619    
620     ContactElement_Nodes[0][k]=Element_Nodes[1][kintold]
621     ContactElement_Nodes[1][k]=Element_Nodes[2][kintold]
622     ContactElement_Nodes[2][k]=Element_Nodes[6][kintold]
623     ContactElement_Nodes[3][k]=Element_Nodes[5][kintold]
624     ContactElement_Nodes[4][k]=Element_Nodes[0][kintold]
625     ContactElement_Nodes[5][k]=Element_Nodes[3][kintold]
626     ContactElement_Nodes[6][k]=Element_Nodes[7][kintold]
627     ContactElement_Nodes[7][k]=Element_Nodes[4][kintold]
628    
629     ContactElement_Nodes[8][k]=Element_Nodes[0][kintfault]
630     ContactElement_Nodes[9][k]=Element_Nodes[3][kintfault]
631     ContactElement_Nodes[10][k]=Element_Nodes[7][kintfault]
632     ContactElement_Nodes[11][k]=Element_Nodes[4][kintfault]
633     ContactElement_Nodes[12][k]=Element_Nodes[1][kintfault]
634     ContactElement_Nodes[13][k]=Element_Nodes[2][kintfault]
635     ContactElement_Nodes[14][k]=Element_Nodes[6][kintfault]
636     ContactElement_Nodes[15][k]=Element_Nodes[5][kintfault]
637    
638     meshfaultL.write("%s %d\n"%(ContactElement_Type,ContactElement_Num))
639    
640     for i in range(ContactElement_Num):
641     meshfaultL.write("%d %d"%(ContactElement_ref[i],ContactElement_tag[i]))
642     for j in range(ContactElement_numNodes):
643     meshfaultL.write(" %d"%ContactElement_Nodes[j][i])
644     meshfaultL.write("\n")
645    
646     # point sources (not supported yet)
647    
648     meshfaultL.write("Point1 0")
649    
650    
651    
652     meshfaultL.close()
653    
654 gross 893
655     ne_w=int((ne/height)*width+0.5)
656     mydomainfile = faultL(width,width, height,ne, ne, ne_w,contact=True,xstart=fstart,xend=fend)

  ViewVC Help
Powered by ViewVC 1.1.26