← Back to branch summary

~vcs-imports/escript-finley/trunk

« back to all changes in this revision

Viewing changes to dudley/test/python/slip_stress_mesh_old.py

  • Committer: jfenwick
  • Date: 2010-10-11 01:48:14 UTC
  • Revision ID: svn-v4:77569008-7704-0410-b7a0-a92fef0b09fd:trunk:3259
Merging dudley and scons updates from branches

Show diffs side-by-side

added added

removed removed

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