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

Contents of /trunk/dudley/test/python/slip_stress_mesh_old.py

Parent Directory Parent Directory | Revision Log Revision Log


Revision 3774 - (show annotations)
Wed Jan 18 06:29:34 2012 UTC (7 years, 9 months ago) by jfenwick
File MIME type: text/x-python
File size: 24695 byte(s)
dudley, pasowrap, pycad

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)

  ViewVC Help
Powered by ViewVC 1.1.26