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

Diff of /trunk/finley/test/python/slip_stress_mesh.py

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

revision 892 by gross, Tue Nov 7 08:31:26 2006 UTC revision 893 by gross, Wed Nov 8 08:20:19 2006 UTC
# Line 1  Line 1 
1  # $Id$  # $Id$
                                                                                                                                                                                                       
2  """  """
3    
4  generates   finley mesh simple vertical fault  generates   finley mesh simple vertical fault
# Line 30  from esys.escript import * Line 29  from esys.escript import *
29  from numarray import zeros,Float,array,size  from numarray import zeros,Float,array,size
30    
31  #... generate domain ...  #... generate domain ...
32  ne = 18  ne = 10
33  width  = 100000.  width  = 100000.
34  height =  30000.  height =  30000.
35  fstart=numarray.array([width/2.,7.*width/16.,3.*height/8.])  fstart=numarray.array([width/2.,7.*width/16.,3.*height/8.])
36  fend=numarray.array([width/2.,9.*width/16.,5.*height/8.])  fend=numarray.array([width/2.,9.*width/16.,5.*height/8.])
37    
38  def faultL(l0,l1, l2,ne0, ne1, ne2,contact=False,xstart=zeros(3),xend=zeros(3)):  def faultL(l0,l1, l2,ne0, ne1, ne2,contact=False,xstart=zeros(3),xend=zeros(3)):
    # if contact = true then there is a fault surface  
    # xstart is the co-ordinates of start of fault  
    # xend is the co-ordinates at end of fault  
39     meshfaultL=open('meshfault3D.fly','w')     meshfaultL=open('meshfault3D.fly','w')
40    
41     FaultError1="ERROR: fault defined on or too close to an outer surface"     FaultError1="ERROR: fault defined on or too close to an outer surface"
# Line 111  def faultL(l0,l1, l2,ne0, ne1, ne2,conta Line 107  def faultL(l0,l1, l2,ne0, ne1, ne2,conta
107        n0double=int(i0end)-int(i0start)        n0double=int(i0end)-int(i0start)
108        n1double=int(i1end)-int(i1start)        n1double=int(i1end)-int(i1start)
109        n2double=int(i2end)-int(i2start)        n2double=int(i2end)-int(i2start)
       print "fstart = ",[i0start*l0/ne0, i1start*l1/ne1, i2start*l2/ne2]  
       print "fend = ", [i0start*l0/ne0, i1start*l1/ne1 + n1double*l1/ne1,i2start*l2/ne2 + n2double*l2/ne2 -height]  
110        if (i0start == 0) or (i1start==0) or (i2start==0):        if (i0start == 0) or (i1start==0) or (i2start==0):
111           raise FaultError1           raise FaultError1
112                    
# Line 121  def faultL(l0,l1, l2,ne0, ne1, ne2,conta Line 115  def faultL(l0,l1, l2,ne0, ne1, ne2,conta
115    
116        if n0double==0:        if n0double==0:
117           numNodes=N0*N1*N2+(n1double-1)*(n2double-1)           numNodes=N0*N1*N2+(n1double-1)*(n2double-1)
118            
119        elif n1double==0:        elif n1double==0:
120           numNodes=N0*N1*N2+(n0double-1)*(n2double-1)           numNodes=N0*N1*N2+(n0double-1)*(n2double-1)
121            
# Line 143  def faultL(l0,l1, l2,ne0, ne1, ne2,conta Line 137  def faultL(l0,l1, l2,ne0, ne1, ne2,conta
137     for i2 in range(N2):     for i2 in range(N2):
138        for i1 in range (N1):        for i1 in range (N1):
139           for i0 in range(N0):           for i0 in range(N0):
140              k=M0*i0+M1*i1+M2*i2;              k=  i0 + N0*i1 + N0*N1*i2 # M0*i0+M1*i1+M2*i2;
141              Node_ref[k]=i0 + N0*i1 + N0*N1*i2              Node_ref[k]= i0 + N0*i1 + N0*N1*i2
142              # no periodic boundary conditions              # no periodic boundary conditions
143              Node_DOF[k]=Node_ref[k]              Node_DOF[k]=Node_ref[k]
144              Node_tag[k]=0              Node_tag[k]=0
145              Node[0][k]=(i0)*l0/(N0-1)              Node[0][k]=(i0)*l0/(N0-1)
146              Node[1][k]=(i1)*l1/(N1-1)              Node[1][k]=(i1)*l1/(N1-1)
147              Node[2][k]=(i2)*l2/(N2-1)-height              Node[2][k]=(i2)*l2/(N2-1)
148    
149     # define double nodes on fault (will have same coordinates as some of nodes already defined)     # define double nodes on fault (will have same coordinates as some of nodes already defined)
150     # only get double nodes on INTERIOR of fault     # only get double nodes on INTERIOR of fault
# Line 167  def faultL(l0,l1, l2,ne0, ne1, ne2,conta Line 161  def faultL(l0,l1, l2,ne0, ne1, ne2,conta
161            
162          for i2 in range(n2double-1):          for i2 in range(n2double-1):
163             for i1 in range(n1double-1):             for i1 in range(n1double-1):
164                k=Fault_NE+M1f*i1+M2f*i2                # CHANGED:
165                Node_ref[k]=Fault_NE + i1 + (n1double-1)*i2                k=Fault_NE+i1+(n1double-1)*i2
166                  Node_ref[k]= k #Fault_NE + i1 + (n1double-1)*i2
167                Node_DOF[k]=Node_ref[k]                Node_DOF[k]=Node_ref[k]
168                Node_tag[k]=1                Node_tag[k]=1
169                Node[0][k]=i0start*l0/ne0                Node[0][k]=i0start*l0/ne0
170                Node[1][k]=i1start*l1/ne1 + (i1+1)*l1/ne1                Node[1][k]=i1start*l1/ne1 + (i1+1)*l1/ne1
171                Node[2][k]=i2start*l2/ne2 + (i2+1)*l2/ne2 -height                Node[2][k]=i2start*l2/ne2 + (i2+1)*l2/ne2
   
172        # elif n1double==0:        # elif n1double==0:
173        # elif n2double==0:        # elif n2double==0:
174          print "fstart = ",[i0start*l0/ne0, i1start*l1/ne1                  , i2start*l2/ne2]
175          print "fend = ", [i0start*l0/ne0 , i1start*l1/ne1 + n1double*l1/ne1, i2start*l2/ne2 + n2double*l2/ne2]
176    
177     # write nodes to file     # write nodes to file
178     for i in range(numNodes):     for i in range(numNodes):
# Line 197  def faultL(l0,l1, l2,ne0, ne1, ne2,conta Line 193  def faultL(l0,l1, l2,ne0, ne1, ne2,conta
193    
194     #print 'Interior elements'     #print 'Interior elements'
195    
196       print "M0,M1,M2",M0,M1,M2
197    
198     for i2 in range(ne2):     for i2 in range(ne2):
199        for i1 in range (ne1):        for i1 in range (ne1):
200           for i0 in range(ne0):           for i0 in range(ne0):
# Line 206  def faultL(l0,l1, l2,ne0, ne1, ne2,conta Line 204  def faultL(l0,l1, l2,ne0, ne1, ne2,conta
204              Element_ref[k]=k              Element_ref[k]=k
205              Element_tag[k]=0              Element_tag[k]=0
206              # for hex8 the interior elements are specified by 8 nodes              # for hex8 the interior elements are specified by 8 nodes
207                #CHANGED:
208              Element_Nodes[0][k]=node0;              Element_Nodes[0][k]=node0;
209              Element_Nodes[1][k]=node0+1;              Element_Nodes[1][k]=node0+1;
210              Element_Nodes[2][k]=node0+N0+1;              Element_Nodes[2][k]=node0+N0+1;
# Line 223  def faultL(l0,l1, l2,ne0, ne1, ne2,conta Line 222  def faultL(l0,l1, l2,ne0, ne1, ne2,conta
222           x0e= i0end*l0/ne0           x0e= i0end*l0/ne0
223           x1e= i1end*l1/ne1           x1e= i1end*l1/ne1
224           x2e= i2end*l2/ne2           x2e= i2end*l2/ne2
225           #print x0s,x1s,x2s,x0e,x1e,x2e           #print "x0s,x1s,x2s,x0e,x1e,x2e",  x0s,x1s,x2s,x0e,x1e,x2e
           
226           if (n1double==1) or (n2double==1):           if (n1double==1) or (n2double==1):
227              raise FaultError2              raise FaultError2
228           for i2 in range(n2double):           for i2 in range(n2double):
# Line 232  def faultL(l0,l1, l2,ne0, ne1, ne2,conta Line 230  def faultL(l0,l1, l2,ne0, ne1, ne2,conta
230                 # here the coordinates of kfault and kold are the same                 # here the coordinates of kfault and kold are the same
231                 # Ref for fault node (only on interior nodes of fault):                 # Ref for fault node (only on interior nodes of fault):
232                 if (i1>0) and (i2>0):                 if (i1>0) and (i2>0):
233                    kfault=Fault_NE+M1f*(i1-1.)+M2f*(i2-1.)                    kfault=Fault_NE+(i1-1.)+(n1double-1)*(i2-1.)
234                    #print kfault, Node[0][kfault],Node[1][kfault],Node[2][kfault]                    #print kfault , Node[0][int(kfault)],Node[1][int(kfault)],Node[2][int(kfault)]
235                 else:                 else:
236                    kfault=0.                    kfault=0.
237                 # determine bottom corner node of each element                 # determine bottom corner node of each element
238    
239                 # Ref for normal interior node:                 # Ref for normal interior node:
240                 kold=int(M0*i0start+M1*(i1start + i1) + M2*(i2start+i2))                 kold=int(i0start+N0*(i1start + i1) + (N0*N1)*(i2start+i2))
241                 #print kold, Node[0][kold],Node[1][kold],Node[2][kold]                 #print kold, Node[0][kold],Node[1][kold],Node[2][kold]
242                 # Ref for interior element:                 # Ref for interior element:
243                 kint=int(M0i*i0start + M1i*(i1start+i1) + M2i*(i2start+i2))                           kint=int(i0start + ne0*(i1start+i1) + (ne0*ne1)*(i2start+i2))          
244                 #print kint, Element_Nodes[0][kint]                 #print kint, Element_Nodes[0][kint]
245    
246                 x0= (i0start)*l0/ne0                 x0= (i0start)*l0/ne0
# Line 256  def faultL(l0,l1, l2,ne0, ne1, ne2,conta Line 254  def faultL(l0,l1, l2,ne0, ne1, ne2,conta
254                 # are on the fault:                 # are on the fault:
255                 if (i1==0) and (i2==0):                 if (i1==0) and (i2==0):
256                    # nearest fault node:                    # nearest fault node:
257                    kfaultref=int(Fault_NE+M1f*i1+M2f*i2)                    kfaultref=int(Fault_NE+i1+(n1double-1)*i2)
258                 elif (i1==0):                 elif (i1==0):
259                    # nearest fault node                    # nearest fault node
260                    kfaultref=int(Fault_NE+i1*M1f+(i2-1.)*(M2f))                    kfaultref=int(Fault_NE+i1+(i2-1.)*(n1double-1))
261                 elif (i2==0):                 elif (i2==0):
262                    # nearest fault node                    # nearest fault node
263                    kfaultref=int(Fault_NE+(i1-1.)*M1f + i2*M2f)                    kfaultref=int(Fault_NE+(i1-1.) + i2*(n1double-1))
264                 else:                 else:
265                    # looking at element with fault node on bottom corner                    # looking at element with fault node on bottom corner
266                    kfaultref=int(kfault)                    kfaultref=int(kfault)
# Line 273  def faultL(l0,l1, l2,ne0, ne1, ne2,conta Line 271  def faultL(l0,l1, l2,ne0, ne1, ne2,conta
271    
272                 # overwrite 4 outer corner elements of fault (only one node changed)                 # overwrite 4 outer corner elements of fault (only one node changed)
273                 if (i1==0 and i2==0):                           if (i1==0 and i2==0):          
274                      #nodecheck=int(Element_Nodes[7][kint] )
275                      #print nodecheck, Node[0][nodecheck],Node[1][nodecheck],Node[2][nodecheck]
276    
277                    Element_Nodes[7][kint]=kfaultref                    Element_Nodes[7][kint]=kfaultref
278    
279                      #print kfaultref, Node[0][kfaultref],Node[1][kfaultref],Node[2][kfaultref]
280    
281    
282                 elif (i1==0 and i2==n2double-1):                 elif (i1==0 and i2==n2double-1):
283    
284                      #nodecheck=int(Element_Nodes[3][kint] )
285                      #print nodecheck, Node[0][nodecheck],Node[1][nodecheck],Node[2][nodecheck]
286    
287                    Element_Nodes[3][kint]=kfaultref                    Element_Nodes[3][kint]=kfaultref
288    
289                      #print kfaultref, Node[0][kfaultref],Node[1][kfaultref],Node[2][kfaultref]
290    
291                 elif (i1==n1double-1 and i2==0):                 elif (i1==n1double-1 and i2==0):
292                      #nodecheck=int(Element_Nodes[4][kint] )
293                      #print nodecheck, Node[0][nodecheck],Node[1][nodecheck],Node[2][nodecheck]
294    
295                    Element_Nodes[4][kint]=kfaultref                    Element_Nodes[4][kint]=kfaultref
296                      #print kfaultref, Node[0][kfaultref],Node[1][kfaultref],Node[2][kfaultref]
297    
298                 elif (i1==n1double-1 and i2==n2double-1):                 elif (i1==n1double-1 and i2==n2double-1):
299                      #nodecheck=int(Element_Nodes[0][kint] )
300                      #print nodecheck, Node[0][nodecheck],Node[1][nodecheck],Node[2][nodecheck]
301    
302                    Element_Nodes[0][kint]=kfaultref                    Element_Nodes[0][kint]=kfaultref
303                      #print kfaultref, Node[0][kfaultref],Node[1][kfaultref],Node[2][kfaultref]
304    
305                 # overwrite 4 sides of fault (only 2 nodes changed)                 # overwrite 4 sides of fault (only 2 nodes changed)
306                 elif (i1==0):                 elif (i1==0):
307    
308                      #nodecheck=int(Element_Nodes[3][kint] )
309                      #print nodecheck, Node[0][nodecheck],Node[1][nodecheck],Node[2][nodecheck]
310                      #nodecheck=int(Element_Nodes[7][kint] )
311                      #print nodecheck, Node[0][nodecheck],Node[1][nodecheck],Node[2][nodecheck]
312    
313                    Element_Nodes[3][kint]=kfaultref                    Element_Nodes[3][kint]=kfaultref
314                    kfaultref1=int(kfaultref+M2f)                    kfaultref1=int(kfaultref+(n1double-1))
315                    Element_Nodes[7][kint]=kfaultref1                    Element_Nodes[7][kint]=kfaultref1
316    
317                      #print kfaultref, Node[0][kfaultref],Node[1][kfaultref],Node[2][kfaultref]
318                      #print kfaultref1, Node[0][kfaultref1],Node[1][kfaultref1],Node[2][kfaultref1]
319    
320    
321    
322                 elif (i1==n1double-1):                 elif (i1==n1double-1):
323    
324                      #nodecheck=int(Element_Nodes[0][kint] )
325                      #print nodecheck, Node[0][nodecheck],Node[1][nodecheck],Node[2][nodecheck]
326                      #nodecheck=int(Element_Nodes[4][kint] )
327                      #print nodecheck, Node[0][nodecheck],Node[1][nodecheck],Node[2][nodecheck]
328    
329    
330                    Element_Nodes[0][kint]=kfaultref                                        Element_Nodes[0][kint]=kfaultref                    
331                    kfaultref1=kfaultref+M2f                    kfaultref1=kfaultref+(n1double-1)
332                    Element_Nodes[4][kint]=kfaultref1                    Element_Nodes[4][kint]=kfaultref1
333    
334                      #print kfaultref, Node[0][kfaultref],Node[1][kfaultref],Node[2][kfaultref]
335                      #print kfaultref1, Node[0][kfaultref1],Node[1][kfaultref1],Node[2][kfaultref1]
336    
337                 elif (i2==0):                 elif (i2==0):
338    
339                      #nodecheck=int(Element_Nodes[4][kint] )
340                      #print nodecheck, Node[0][nodecheck],Node[1][nodecheck],Node[2][nodecheck]
341                      #nodecheck=int(Element_Nodes[7][kint] )
342                      #print nodecheck, Node[0][nodecheck],Node[1][nodecheck],Node[2][nodecheck]
343    
344                    Element_Nodes[4][kint]=kfaultref                    Element_Nodes[4][kint]=kfaultref
345                    kfaultref1=kfaultref+M1f                    kfaultref1=kfaultref+1
346                    Element_Nodes[7][kint]=kfaultref1                    Element_Nodes[7][kint]=kfaultref1
347    
                elif (i2==n2double-1):  
                   #print i1  
                   #nodecheck=int(Element_Nodes[0][kint] )  
                   #print nodecheck, Node[0][nodecheck],Node[1][nodecheck],Node[2][nodecheck]  
348                    #print kfaultref, Node[0][kfaultref],Node[1][kfaultref],Node[2][kfaultref]                    #print kfaultref, Node[0][kfaultref],Node[1][kfaultref],Node[2][kfaultref]
349                      #print kfaultref1, Node[0][kfaultref1],Node[1][kfaultref1],Node[2][kfaultref1]
350    
351                    Element_Nodes[0][kint]=kfaultref                 elif (i2==n2double-1):
                   kfaultref1=kfaultref+M1f  
352    
353                      #nodecheck=int(Element_Nodes[0][kint] )
354                      #print nodecheck, Node[0][nodecheck],Node[1][nodecheck],Node[2][nodecheck]
355                    #nodecheck=int(Element_Nodes[3][kint] )                    #nodecheck=int(Element_Nodes[3][kint] )
356                    #print nodecheck, Node[0][nodecheck],Node[1][nodecheck],Node[2][nodecheck]                    #print nodecheck, Node[0][nodecheck],Node[1][nodecheck],Node[2][nodecheck]
                   #print kfaultref1, Node[0][kfaultref1],Node[1][kfaultref1],Node[2][kfaultref1]  
357    
358                      Element_Nodes[0][kint]=kfaultref
359                      kfaultref1=kfaultref+1
360                    Element_Nodes[3][kint]=kfaultref1                    Element_Nodes[3][kint]=kfaultref1
361    
362                      #print kfaultref, Node[0][kfaultref],Node[1][kfaultref],Node[2][kfaultref]
363                      #print kfaultref1, Node[0][kfaultref1],Node[1][kfaultref1],Node[2][kfaultref1]
364    
365                 # overwrite interior fault elements (4 nodes changed)                 # overwrite interior fault elements (4 nodes changed)
366                 else:                 else:
                   #print i1,i2  
367                    #nodecheck=int(Element_Nodes[0][kint] )                    #nodecheck=int(Element_Nodes[0][kint] )
368                    #print nodecheck, Node[0][nodecheck],Node[1][nodecheck],Node[2][nodecheck]                    #print nodecheck, Node[0][nodecheck],Node[1][nodecheck],Node[2][nodecheck]
369                    #print kfaultref, Node[0][kfaultref],Node[1][kfaultref],Node[2][kfaultref]                    #print kfaultref, Node[0][kfaultref],Node[1][kfaultref],Node[2][kfaultref]
370    
371                    Element_Nodes[0][kint]=kfaultref                    Element_Nodes[0][kint]=kfaultref
372                    #if (x1<x1e and x2<x2e):                    #if (x1<x1e and x2<x2e):
373                    kfaultref1=kfaultref+M1f                    kfaultref1=kfaultref+1
374                    kfaultref2=kfaultref+M2f                    kfaultref2=kfaultref+(n1double-1)
375                    kfaultref3=kfaultref+M1f+M2f                    kfaultref3=kfaultref+1+(n1double-1)
376    
377                    nodecheck=int(Element_Nodes[3][kint] )                    #nodecheck=int(Element_Nodes[3][kint] )
378                    #print nodecheck, Node[0][nodecheck],Node[1][nodecheck],Node[2][nodecheck]                    #print nodecheck, Node[0][nodecheck],Node[1][nodecheck],Node[2][nodecheck]
379                    #print kfaultref1, Node[0][kfaultref1],Node[1][kfaultref1],Node[2][kfaultref1]                    #print kfaultref1, Node[0][kfaultref1],Node[1][kfaultref1],Node[2][kfaultref1]
380    
381                    nodecheck=int(Element_Nodes[4][kint] )                    #nodecheck=int(Element_Nodes[4][kint] )
382                    #print nodecheck, Node[0][nodecheck],Node[1][nodecheck],Node[2][nodecheck]                    #print nodecheck, Node[0][nodecheck],Node[1][nodecheck],Node[2][nodecheck]
383                    #print kfaultref2, Node[0][kfaultref2],Node[1][kfaultref2],Node[2][kfaultref2]                    #print kfaultref2, Node[0][kfaultref2],Node[1][kfaultref2],Node[2][kfaultref2]
384    
385                    nodecheck=int(Element_Nodes[7][kint] )                    #nodecheck=int(Element_Nodes[7][kint] )
386                    #print nodecheck, Node[0][nodecheck],Node[1][nodecheck],Node[2][nodecheck]                    #print nodecheck, Node[0][nodecheck],Node[1][nodecheck],Node[2][nodecheck]
387                    #print kfaultref3, Node[0][kfaultref3],Node[1][kfaultref3],Node[2][kfaultref3]                    #print kfaultref3, Node[0][kfaultref3],Node[1][kfaultref3],Node[2][kfaultref3]
388    
# Line 527  def faultL(l0,l1, l2,ne0, ne1, ne2,conta Line 572  def faultL(l0,l1, l2,ne0, ne1, ne2,conta
572                 #print k                 #print k
573                 # define reference for interior elements with x0<=x0s                 # define reference for interior elements with x0<=x0s
574                 # here the nodes are the old interior nodes                 # here the nodes are the old interior nodes
575                 kintold=int(M0i*(i0start-1) + M1i*(i1start+i1) + M2i*(i2start+i2))                 kintold=int((i0start-1) + ne0*(i1start+i1) + ne0*ne1*(i2start+i2))
576    
577                 # define reference for interior elements with x0>x0s                 # define reference for interior elements with x0>x0s
578                 # here the double nodes are the fault nodes                 # here the double nodes are the fault nodes
579    
580                 kintfault=int(M0i*i0start + M1i*(i1start+i1) + M2i*(i2start+i2))                 kintfault=int(i0start + ne0*(i1start+i1) + ne0*ne1*(i2start+i2))
581    
582                 #nodecheck=int(Element_Nodes[1][kintold] )                 #nodecheck=int(Element_Nodes[1][kintold] )
583                 #print nodecheck, Node[0][nodecheck],Node[1][nodecheck],Node[2][nodecheck]                 #print nodecheck, Node[0][nodecheck],Node[1][nodecheck],Node[2][nodecheck]
# Line 591  def faultL(l0,l1, l2,ne0, ne1, ne2,conta Line 636  def faultL(l0,l1, l2,ne0, ne1, ne2,conta
636    
637     meshfaultL.close()     meshfaultL.close()
638    
639  # ne_w=int((ne/height)*width+0.5)  
640  ne_w=ne  ne_w=int((ne/height)*width+0.5)
641  mydomainfile = faultL(width,width, height,ne, ne, ne,contact=True,xstart=fstart,xend=fend)  mydomainfile = faultL(width,width, height,ne, ne, ne_w,contact=True,xstart=fstart,xend=fend)
Legend:
Removed from v.892  
changed lines
  Added in v.893
  ViewVC Help
Powered by ViewVC 1.1.26