/[escript]/trunk/doc/examples/cookbook/example09n.py
ViewVC logotype

Diff of /trunk/doc/examples/cookbook/example09n.py

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

revision 3388 by jfenwick, Mon Oct 11 01:48:14 2010 UTC revision 3389 by ahallam, Wed Dec 1 23:03:35 2010 UTC
# Line 35  from esys.pycad.gmsh import Design #Fini Line 35  from esys.pycad.gmsh import Design #Fini
35  from esys.finley import MakeDomain #Converter for escript  from esys.finley import MakeDomain #Converter for escript
36  from esys.escript import mkDir, getMPISizeWorld  from esys.escript import mkDir, getMPISizeWorld
37  import os  import os
38    import math
39  ########################################################MPI WORLD CHECK  ########################################################MPI WORLD CHECK
40  if getMPISizeWorld() > 1:  if getMPISizeWorld() > 1:
41      import sys      import sys
# Line 42  if getMPISizeWorld() > 1: Line 43  if getMPISizeWorld() > 1:
43      sys.exit(0)      sys.exit(0)
44    
45  # make sure path exists  # make sure path exists
46  save_path= os.path.join("data","example09")  save_path= os.path.join("data","example09c")
47  mkDir(save_path)  mkDir(save_path)
48    
49  ################################################ESTABLISHING PARAMETERS  ################################################ESTABLISHING PARAMETERS
50  #Model Parameters  #Model Parameters
51  xwidth=2000.0   #x width of model  origin=[0,0]                  #orign of model
52  ywidth=2000.0   #y width of model  xwidth=300.0                      #width of model
53  depth=500.0   #depth of model  depth=-100.0                      #depth of model
54  intf=depth/2.   #Depth of the interface.  nintf=3                         #number of the interfaces
55    lintf_depths=[-20,-40,-60]       #depth of interfaces
56    rintf_depths=[-30,-50,-70]      #vertical displacement across fault
57    fault_dip=40.0                  #dip of fault plane
58    fault_atsurface=50.0             #location of fault at surface
59    fault_width=10                 #apparent width of fault plane
60    
61    element_size=1.
62    
63  ####################################################DOMAIN CONSTRUCTION  ####################################################DOMAIN CONSTRUCTION
64  # Domain Corners  x0=0.0+origin[0]
65  p0=Point(0.0,    0.0,      0.0)  y0=0.0+origin[1]
66  p1=Point(xwidth, 0.0,      0.0)  #z=0.0+origin[2]
67  p2=Point(xwidth, ywidth,   0.0)  fault_atsurface=fault_atsurface+origin[0]
68  p3=Point(0.0,    ywidth,   0.0)  xwidth=xwidth+origin[0]
69  p4=Point(0.0,    ywidth, depth)  depth=depth+origin[1]
70  p5=Point(0.0,    0.0,    depth)  # Construction points, 4 vectors that descend from the surface with nintf+2 points.
71  p6=Point(xwidth, 0.0,    depth)  left_edge=[Point(x0,y0)];
72  p7=Point(xwidth, ywidth, depth)  leftf_edge=[Point(fault_atsurface,y0)];
73  # Join corners in anti-clockwise manner.  rightf_edge=[Point(fault_atsurface+fault_width,y0)];
74  l01=Line(p0, p1)  right_edge=[Point(xwidth,y0)];
75  l12=Line(p1, p2)  
76  l23=Line(p2, p3)  for i in range(0,nintf):
77  l30=Line(p3, p0)      left_shift=(lintf_depths[i]-y0)/math.tan(fault_dip)
78  l56=Line(p5, p6)      right_shift=(rintf_depths[i]-y0)/math.tan(fault_dip)
79  l67=Line(p6, p7)      left_edge.append(Point(x0,lintf_depths[i]+origin[1]))
80  l74=Line(p7, p4)      leftf_edge.append(Point(fault_atsurface+left_shift,lintf_depths[i]+origin[1]))
81  l45=Line(p4, p5)      rightf_edge.append(Point(fault_atsurface+right_shift+fault_width,rintf_depths[i]+origin[1]))
82        right_edge.append(Point(xwidth,rintf_depths[i]+origin[1]))
83  # Join line segments to create domain boundaries and then surfaces  
84  ctop=CurveLoop(l01, l12, l23, l30);     stop=PlaneSurface(ctop)  left_edge.append(Point(x0,depth))
85  cbot=CurveLoop(-l67, -l56, -l45, -l74); sbot=PlaneSurface(cbot)  leftf_edge.append(Point(fault_atsurface+(depth-y0)/math.tan(fault_dip),depth))
86    rightf_edge.append(Point(fault_atsurface+fault_width+(depth-y0)/math.tan(fault_dip),depth))
87  # for each side  right_edge.append(Point(xwidth,depth))
88  ip0=Point(0.0,    0.0,      intf)  
89  ip1=Point(xwidth, 0.0,      intf)  #Build lines
90  ip2=Point(xwidth, ywidth,   intf)  lright=[]; nlright=[];
91  ip3=Point(0.0,    ywidth,   intf)  lhright=[]; nlhright=[];
92    lfright=[]; nlfright=[];
93  linte_ar=[]; #lines for vertical edges  lfhor=[]; nlfhor=[];
94  linhe_ar=[]; #lines for horizontal edges  lfleft=[]; nlfleft=[];
95  linte_ar.append(Line(p0,ip0))  lhleft=[]; nlhleft=[];
96  linte_ar.append(Line(ip0,p5))  lleft=[]; nlleft=[];
97  linte_ar.append(Line(p1,ip1))  
98  linte_ar.append(Line(ip1,p6))  #Build vertical lines
99  linte_ar.append(Line(p2,ip2))  for i in range(0,nintf+1):
100  linte_ar.append(Line(ip2,p7))      lleft.append(Line(left_edge[i],left_edge[i+1]))
101  linte_ar.append(Line(p3,ip3))      lfleft.append(Line(leftf_edge[i],leftf_edge[i+1]))
102  linte_ar.append(Line(ip3,p4))      lfright.append(Line(rightf_edge[i],rightf_edge[i+1]))
103        lright.append(Line(right_edge[i],right_edge[i+1]))
104  linhe_ar.append(Line(ip0,ip1))  
105  linhe_ar.append(Line(ip1,ip2))  #Build horizontal lines
106  linhe_ar.append(Line(ip2,ip3))  for i in range(0,nintf+2):
107  linhe_ar.append(Line(ip3,ip0))      lhleft.append(Line(left_edge[i],leftf_edge[i]))
108        lhright.append(Line(rightf_edge[i],right_edge[i]))
109  cintfa_ar=[]; cintfb_ar=[] #curveloops for above and below interface on sides  lfhor.append(Line(leftf_edge[0],rightf_edge[0]))
110  cintfa_ar.append(CurveLoop(linte_ar[0],linhe_ar[0],-linte_ar[2],-l01))  lfhor.append(Line(leftf_edge[nintf+1],rightf_edge[nintf+1]))
111  cintfa_ar.append(CurveLoop(linte_ar[2],linhe_ar[1],-linte_ar[4],-l12))  
112  cintfa_ar.append(CurveLoop(linte_ar[4],linhe_ar[2],-linte_ar[6],-l23))  #Build negative lines
113  cintfa_ar.append(CurveLoop(linte_ar[6],linhe_ar[3],-linte_ar[0],-l30))  for i in range(nintf,-1,-1):
114        nlleft.append(-lleft[i])
115  cintfb_ar.append(CurveLoop(linte_ar[1],l56,-linte_ar[3],-linhe_ar[0]))      nlfleft.append(-lfleft[i])
116  cintfb_ar.append(CurveLoop(linte_ar[3],l67,-linte_ar[5],-linhe_ar[1]))      nlfright.append(-lfright[i])
117  cintfb_ar.append(CurveLoop(linte_ar[5],l74,-linte_ar[7],-linhe_ar[2]))      nlright.append(-lright[i])
118  cintfb_ar.append(CurveLoop(linte_ar[7],l45,-linte_ar[1],-linhe_ar[3]))  for i in range(nintf+1,-1,-1):
119        nlhleft.append(-lhleft[i])
120  sintfa_ar=[PlaneSurface(cintfa_ar[i]) for i in range(0,4)]      nlhright.append(-lhright[i])
121  sintfb_ar=[PlaneSurface(cintfb_ar[i]) for i in range(0,4)]  
122    #Build curveloops
123  sintf=PlaneSurface(CurveLoop(*tuple(linhe_ar)))  lcurves=[]
124    fcurves=[]
125  vintfa=Volume(SurfaceLoop(stop,-sintf,*tuple(sintfa_ar)))  rcurves=[]
126  vintfb=Volume(SurfaceLoop(sbot,sintf,*tuple(sintfb_ar)))  
127    #Fault
128  # Create the volume.  for i in range(0,nintf+1):
129  #sloop=SurfaceLoop(stop,sbot,*tuple(sintfa_ar+sintfb_ar))      fcurves.append(lfleft[i])
130  #model=Volume(sloop)  fcurves.append(lfhor[1])
131    for i in range(0,nintf+1):
132        fcurves.append(nlfright[i])
133  #############################################EXPORTING MESH FOR ESCRIPT  fcurves.append(-lfhor[0])
134  # Create a Design which can make the mesh  fcurves=CurveLoop(*tuple(fcurves))
135  d=Design(dim=3, element_size=200.0)  
136    #Left and Right Blocks
137  d.addItems(PropertySet('vintfa',vintfa))  for i in range(0,nintf+1):
138  d.addItems(PropertySet('vintfb',vintfb))      lcurves.append(CurveLoop(lleft[i],lhleft[i+1],nlfleft[nintf-i],nlhleft[nintf+1-i]))
139  d.addItems(sintf)      rcurves.append(CurveLoop(lfright[i],lhright[i+1],nlright[nintf-i],nlhright[nintf+1-i]))
140        
141    #Build Surfaces
142    fsurf=PlaneSurface(fcurves)
143    lsurf=[]
144    rsurf=[]
145    for i in range(0,nintf+1):
146        lsurf.append(PlaneSurface(lcurves[i]))
147        rsurf.append(PlaneSurface(rcurves[i]))
148    
149    d=Design(dim=2, element_size=element_size, order=2)
150    
151    d.addItems(PropertySet('fault',fsurf))
152    for i in range(0,nintf+1):
153        d.addItems(PropertySet('lblock%d'%i,lsurf[i]))
154        d.addItems(PropertySet('rblock%d'%i,rsurf[i]))
155    
156    d.addItems(PropertySet('top',lhright[0],lfhor[0],lhleft[0],lhright[4],lhleft[4],lfhor[1]))
157    #for i in range(0,2):
158    #     d.addItems(lfhor[i])
159    
160    
161    
162  d.setScriptFileName(os.path.join(save_path,"example09n.geo"))  d.setScriptFileName(os.path.join(save_path,"example09n.geo"))
163  d.setMeshFileName(os.path.join(save_path,"example09n.msh"))  d.setMeshFileName(os.path.join(save_path,"example09n.msh"))
164  #  #
165  #    #  make the domain:
166  #  #
167  domain=MakeDomain(d)  domain=MakeDomain(d)
 # Create a file that can be read back in to python with  
168  # mesh=ReadMesh(fileName)  # mesh=ReadMesh(fileName)
169  domain.write(os.path.join(save_path,"example09n.fly"))  domain.write(os.path.join(save_path,"example09n.fly"))
170    
 from esys.pycad import layer_cake  
 intfaces=[10,30,50,55,80,100,200,250,400,500]  
   
 cmplx_domain=layer_cake.LayerCake(xwidth,ywidth,intfaces,200.0)  
 cmplx_domain.setScriptFileName(os.path.join(save_path,"example09lc.geo"))  
 cmplx_domain.setMeshFileName(os.path.join(save_path,"example09lc.msh"))  
 dcmplx=MakeDomain(cmplx_domain)  
 dcmplx.write(os.path.join(save_path,"example09lc.fly"))  
171    
172    
173    

Legend:
Removed from v.3388  
changed lines
  Added in v.3389

  ViewVC Help
Powered by ViewVC 1.1.26