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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 6651 - (show annotations)
Wed Feb 7 02:12:08 2018 UTC (20 months, 1 week ago) by jfenwick
File MIME type: text/x-python
File size: 24862 byte(s)
Make everyone sad by touching all the files

Copyright dates update

1
2 ##############################################################################
3 #
4 # Copyright (c) 2003-2018 by The University of Queensland
5 # http://www.uq.edu.au
6 #
7 # Primary Business: Queensland, Australia
8 # Licensed under the Apache License, version 2.0
9 # http://www.apache.org/licenses/LICENSE-2.0
10 #
11 # Development until 2012 by Earth Systems Science Computational Center (ESSCC)
12 # Development 2012-2013 by School of Earth Sciences
13 # Development from 2014 by Centre for Geoscience Computing (GeoComp)
14 #
15 ##############################################################################
16
17 from __future__ import print_function, division
18
19 __copyright__="""Copyright (c) 2003-2018 by The University of Queensland
20 http://www.uq.edu.au
21 Primary Business: Queensland, Australia"""
22 __license__="""Licensed under the Apache License, version 2.0
23 http://www.apache.org/licenses/LICENSE-2.0"""
24 __url__="https://launchpad.net/escript-finley"
25
26 """
27
28 generates finley mesh simple vertical fault
29
30 THIS CODE CREATES RICH CONTACT ELEMENTS AND RICH FACE ELEMENTS
31 with fix for contact elements at FAULT ENDS
32
33 :var __author__: name of author
34 :var __copyright__: copyrights
35 :var __license__: licence agreement
36 :var __url__: url entry point on documentation
37 :var __version__: version
38 :var __date__: date of the version
39 """
40
41 __author__="Louise Kettle"
42
43 from esys.escript import *
44 from numpy import zeros,float64,array,size
45
46 #... generate domain ...
47 ne = 10
48 width = 100000.
49 height = 30000.
50 fstart=array([width/2.,7.*width/16.,3.*height/8.])
51 fend=array([width/2.,9.*width/16.,5.*height/8.])
52
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
134 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],float64)
144 Node_ref=zeros(numNodes,float64)
145 Node_DOF=zeros(numNodes,float64)
146 Node_tag=zeros(numNodes,float64)
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 k= i0 + N0*i1 + N0*N1*i2 # M0*i0+M1*i1+M2*i2;
156 Node_ref[k]= i0 + N0*i1 + N0*N1*i2
157 # 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 Node[2][k]=(i2)*l2/(N2-1)
163
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 # CHANGED:
180 k=Fault_NE+i1+(n1double-1)*i2
181 Node_ref[k]= k #Fault_NE + i1 + (n1double-1)*i2
182 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 Node[2][k]=i2start*l2/ne2 + (i2+1)*l2/ne2
187 # elif n1double==0:
188 # elif n2double==0:
189 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
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],float64)
206 Element_ref=zeros(ne0*ne1*ne2,float64)
207 Element_tag=zeros(ne0*ne1*ne2,float64)
208
209 #print 'Interior elements'
210
211 print("M0,M1,M2",M0,M1,M2)
212
213 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 #CHANGED:
223 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 #print "x0s,x1s,x2s,x0e,x1e,x2e", x0s,x1s,x2s,x0e,x1e,x2e
241 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 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 else:
251 kfault=0.
252 # determine bottom corner node of each element
253
254 # Ref for normal interior node:
255 kold=int(i0start+N0*(i1start + i1) + (N0*N1)*(i2start+i2))
256 #print kold, Node[0][kold],Node[1][kold],Node[2][kold]
257 # Ref for interior element:
258 kint=int(i0start + ne0*(i1start+i1) + (ne0*ne1)*(i2start+i2))
259 #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 kfaultref=int(Fault_NE+i1+(n1double-1)*i2)
273 elif (i1==0):
274 # nearest fault node
275 kfaultref=int(Fault_NE+i1+(i2-1.)*(n1double-1))
276 elif (i2==0):
277 # nearest fault node
278 kfaultref=int(Fault_NE+(i1-1.) + i2*(n1double-1))
279 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 #nodecheck=int(Element_Nodes[7][kint] )
290 #print nodecheck, Node[0][nodecheck],Node[1][nodecheck],Node[2][nodecheck]
291
292 Element_Nodes[7][kint]=kfaultref
293
294 #print kfaultref, Node[0][kfaultref],Node[1][kfaultref],Node[2][kfaultref]
295
296
297 elif (i1==0 and i2==n2double-1):
298
299 #nodecheck=int(Element_Nodes[3][kint] )
300 #print nodecheck, Node[0][nodecheck],Node[1][nodecheck],Node[2][nodecheck]
301
302 Element_Nodes[3][kint]=kfaultref
303
304 #print kfaultref, Node[0][kfaultref],Node[1][kfaultref],Node[2][kfaultref]
305
306 elif (i1==n1double-1 and i2==0):
307 #nodecheck=int(Element_Nodes[4][kint] )
308 #print nodecheck, Node[0][nodecheck],Node[1][nodecheck],Node[2][nodecheck]
309
310 Element_Nodes[4][kint]=kfaultref
311 #print kfaultref, Node[0][kfaultref],Node[1][kfaultref],Node[2][kfaultref]
312
313 elif (i1==n1double-1 and i2==n2double-1):
314 #nodecheck=int(Element_Nodes[0][kint] )
315 #print nodecheck, Node[0][nodecheck],Node[1][nodecheck],Node[2][nodecheck]
316
317 Element_Nodes[0][kint]=kfaultref
318 #print kfaultref, Node[0][kfaultref],Node[1][kfaultref],Node[2][kfaultref]
319
320 # overwrite 4 sides of fault (only 2 nodes changed)
321 elif (i1==0):
322
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 Element_Nodes[3][kint]=kfaultref
329 kfaultref1=int(kfaultref+(n1double-1))
330 Element_Nodes[7][kint]=kfaultref1
331
332 #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 elif (i1==n1double-1):
338
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 Element_Nodes[0][kint]=kfaultref
346 kfaultref1=kfaultref+(n1double-1)
347 Element_Nodes[4][kint]=kfaultref1
348
349 #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 elif (i2==0):
353
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 Element_Nodes[4][kint]=kfaultref
360 kfaultref1=kfaultref+1
361 Element_Nodes[7][kint]=kfaultref1
362
363 #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 elif (i2==n2double-1):
367
368 #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 Element_Nodes[0][kint]=kfaultref
374 kfaultref1=kfaultref+1
375 Element_Nodes[3][kint]=kfaultref1
376
377 #print kfaultref, Node[0][kfaultref],Node[1][kfaultref],Node[2][kfaultref]
378 #print kfaultref1, Node[0][kfaultref1],Node[1][kfaultref1],Node[2][kfaultref1]
379
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 kfaultref1=kfaultref+1
389 kfaultref2=kfaultref+(n1double-1)
390 kfaultref3=kfaultref+1+(n1double-1)
391
392 #nodecheck=int(Element_Nodes[3][kint] )
393 #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 #nodecheck=int(Element_Nodes[4][kint] )
397 #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 #nodecheck=int(Element_Nodes[7][kint] )
401 #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],float64)
436 FaceElement_ref=zeros(FaceElement_Num,float64)
437 FaceElement_tag=zeros(FaceElement_Num,float64)
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],float64)
580 ContactElement_ref=zeros(ContactElement_Num,float64)
581 ContactElement_tag=zeros(ContactElement_Num,float64)
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 kintold=int((i0start-1) + ne0*(i1start+i1) + ne0*ne1*(i2start+i2))
591
592 # define reference for interior elements with x0>x0s
593 # here the double nodes are the fault nodes
594
595 kintfault=int(i0start + ne0*(i1start+i1) + ne0*ne1*(i2start+i2))
596
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
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