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

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

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

revision 3068 by ahallam, Tue Jul 6 02:10:59 2010 UTC revision 3069 by ahallam, Wed Jul 21 03:24:48 2010 UTC
# Line 26  Author: Antony Hallam antony.hallam@uqco Line 26  Author: Antony Hallam antony.hallam@uqco
26    
27  ############################################################FILE HEADER  ############################################################FILE HEADER
28  # example09m.py  # example09m.py
29  # Create a simple 3D model for use in example09.  # Create a simple 3D model for use in example09. This is the low res
30    # mesh for illustration purposes only.
31  #  #
32  #######################################################EXTERNAL MODULES  #######################################################EXTERNAL MODULES
33  from esys.pycad import * #domain constructor  from esys.pycad import * #domain constructor
34  from esys.pycad.gmsh import Design #Finite Element meshing package  from esys.pycad.gmsh import Design #Finite Element meshing package
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
 from esys.escript.unitsSI import *  
37  import os  import os
 from math import *  
 import pylab as pl  
 import numpy as np  
38  ########################################################MPI WORLD CHECK  ########################################################MPI WORLD CHECK
39  if getMPISizeWorld() > 1:  if getMPISizeWorld() > 1:
40      import sys      import sys
# Line 50  mkDir(save_path) Line 47  mkDir(save_path)
47    
48  ################################################ESTABLISHING PARAMETERS  ################################################ESTABLISHING PARAMETERS
49  #Model Parameters  #Model Parameters
50  xwidth=2000.0*m   #x width of model  xwidth=1000.0   #x width of model
51  ywidth=2000.0*m   #y width of model  ywidth=1000.0   #y width of model
52  depth=500.0*m   #depth of model  depth=200.0   #depth of model
53    intf=depth/2.   #Depth of the interface.
54    element_size=5.0
55    
56  ####################################################DOMAIN CONSTRUCTION  ####################################################DOMAIN CONSTRUCTION
57  # Domain Corners  # Domain Corners
# Line 78  l45=Line(p4, p5) Line 77  l45=Line(p4, p5)
77  ctop=CurveLoop(l01, l12, l23, l30);     stop=PlaneSurface(ctop)  ctop=CurveLoop(l01, l12, l23, l30);     stop=PlaneSurface(ctop)
78  cbot=CurveLoop(-l67, -l56, -l45, -l74); sbot=PlaneSurface(cbot)  cbot=CurveLoop(-l67, -l56, -l45, -l74); sbot=PlaneSurface(cbot)
79    
 #create interface  
 sspl=51 #number of discrete points in spline  
 xdsp=xwidth/(sspl-1) #dx of spline steps for xwidth  
 ydsp=ywidth/(sspl-1) #dy of spline steps for ywidth  
 dep_sp=250.0*m #avg depth of spline  
 h_sp=25.0*m #heigh of spline  
 orit=-1.0 #orientation of spline 1.0=>up -1.0=>down  
   
 # Generate Material Boundary  
 #x=[ Point(i*xdsp\  
  #   ,-dep_sp+orit*h_sp*cos(pi*i*xdsp/dep_sp+pi))\  
   #   for i in range(0,sspl)\  
   #]  
 #function along x  
 x=[1+orit*cos(2*sspl*pi*i*xdsp/xwidth)\  
      for i in range(0,sspl)]  
 #function along y  
 y=[1+orit*cos(2*sspl*pi*i*ydsp/ywidth)\  
      for i in range(0,sspl)]  
 #containers for surface data  
 xs=np.zeros(sspl,'float')  
 ys=np.zeros(sspl,'float')  
 zs=np.zeros([sspl,sspl],'float')  
   
 for i in range(0,sspl):  
     xs[i]=i*xdsp  
     for j in range(0,sspl):  
         if (i == 0):  
             ys[j]=j*ydsp  
         zs[i,j]=x[i]*y[j]*h_sp+dep_sp  
   
 pl.plot(x); pl.plot(y,'ro')  
 pl.savefig(os.path.join(save_path,"interface.png"))  
 pl.clf()  
 pl.contourf(xs,ys,zs)  
 pl.savefig(os.path.join(save_path,"interface_cf.png"))  
   
 spl_ar=[] #interface spline array  
 linex0_ar=[] #interface line array along x  
 linexy_ar=[] #interface line array along x,ymax  
 sintf_ar=[]  #interface surface array  
 nsintf_ar=[]  
 #loop through x  
 for i in range(0,sspl):  
     #create all the points required  
     point_ar=[Point(xs[i],ys[j],zs[i,j]) for j in range(0,sspl)]  
   
     if (i == 0):  
         ip0=point_ar[0]  
         ip3=point_ar[-1]  
     if (i == sspl-1):  
         ip1=point_ar[0]  
         ip2=point_ar[-1]  
   
     #create the spline and store it  
     spl_ar.append(Spline(*tuple(point_ar)))  
     #create the lines along the x axis and x axis at ymax  
     if (i > 0):  
         #print i,xs[i],ys[0],zs[i,0]  
         tp2=point_ar[0]; tp3=point_ar[-1]  
         tl0=Line(cp0,tp2); tl1=Line(cp1,tp3)  
         linex0_ar.append(tl0); linexy_ar.append(tl1)  
     cp0=point_ar[0]; cp1=point_ar[-1]  
           
 for i in range(0,sspl):  
     #create the mini surfaces via curveloops and then ruledsurfaces  
     if (i < sspl-1):  
         #print 'i',i  
         tc0=CurveLoop(linex0_ar[i],spl_ar[i+1],-linexy_ar[i],-spl_ar[i])  
         ntc0=CurveLoop(spl_ar[i],linexy_ar[i],-spl_ar[i+1],-linex0_ar[i])  
         ts0=RuledSurface(tc0); sintf_ar.append(ts0)  
         nts0=RuledSurface(ntc0); nsintf_ar.append(nts0)  
 #create the interface using the surface loop constructor          
 sintf=SurfaceLoop(*tuple(sintf_ar))  
   
80  # for each side  # for each side
81  #ip0=Point(xs[0],ys[0],zs[0,0])  ip0=Point(0.0,    0.0,      intf)
82  #ip1=Point(xs[-1],ys[0],zs[-1,0])  ip1=Point(xwidth, 0.0,      intf)
83  #ip2=Point(xs[-1],ys[-1],zs[-1,-1])  ip2=Point(xwidth, ywidth,   intf)
84  #ip3=Point(xs[0],ys[-1],zs[0,-1])  ip3=Point(0.0,    ywidth,   intf)
85    
86  linte_ar=[]; #lines for vertical edges  linte_ar=[]; #lines for vertical edges
87  linhe_ar=[]; #lines for horizontal edges  linhe_ar=[]; #lines for horizontal edges
# Line 176  linhe_ar.append(Line(ip2,ip3)) Line 100  linhe_ar.append(Line(ip2,ip3))
100  linhe_ar.append(Line(ip3,ip0))  linhe_ar.append(Line(ip3,ip0))
101    
102  cintfa_ar=[]; cintfb_ar=[] #curveloops for above and below interface on sides  cintfa_ar=[]; cintfb_ar=[] #curveloops for above and below interface on sides
103  cintfa_ar.append(CurveLoop(-linte_ar[2],-l01,linte_ar[0],*tuple(linex0_ar)))  cintfa_ar.append(CurveLoop(linte_ar[0],linhe_ar[0],-linte_ar[2],-l01))
104  cintfa_ar.append(CurveLoop(linte_ar[2],spl_ar[-1],-linte_ar[4],-l12))  cintfa_ar.append(CurveLoop(linte_ar[2],linhe_ar[1],-linte_ar[4],-l12))
105  cintfa_ar.append(CurveLoop(-linte_ar[6],-l23,linte_ar[4],*tuple(linexy_ar)))  cintfa_ar.append(CurveLoop(linte_ar[4],linhe_ar[2],-linte_ar[6],-l23))
106  cintfa_ar.append(CurveLoop(linte_ar[6],-spl_ar[0],-linte_ar[0],-l30))  cintfa_ar.append(CurveLoop(linte_ar[6],linhe_ar[3],-linte_ar[0],-l30))
107    
108  #cintfb_ar.append(CurveLoop(linte_ar[1],l56,-linte_ar[3],-linhe_ar[0]))  cintfb_ar.append(CurveLoop(linte_ar[1],l56,-linte_ar[3],-linhe_ar[0]))
109  #cintfb_ar.append(CurveLoop(linte_ar[3],l67,-linte_ar[5],-linhe_ar[1]))  cintfb_ar.append(CurveLoop(linte_ar[3],l67,-linte_ar[5],-linhe_ar[1]))
110  #cintfb_ar.append(CurveLoop(linte_ar[5],l74,-linte_ar[7],-linhe_ar[2]))  cintfb_ar.append(CurveLoop(linte_ar[5],l74,-linte_ar[7],-linhe_ar[2]))
111  #cintfb_ar.append(CurveLoop(linte_ar[7],l45,-linte_ar[1],-linhe_ar[3]))  cintfb_ar.append(CurveLoop(linte_ar[7],l45,-linte_ar[1],-linhe_ar[3]))
   
 cintfb_ar.append(-CurveLoop(linte_ar[3],-l56,-linte_ar[1],*tuple(linex0_ar)))  
 cintfb_ar.append(-CurveLoop(linte_ar[5],-l67,-linte_ar[3],-spl_ar[-1]))  
 cintfb_ar.append(CurveLoop(linte_ar[5],l74,-linte_ar[7],*tuple(linexy_ar)))  
 cintfb_ar.append(CurveLoop(linte_ar[7],l45,-linte_ar[1],spl_ar[0]))  
   
   
112    
113  sintfa_ar=[PlaneSurface(cintfa_ar[i]) for i in range(0,4)]  sintfa_ar=[PlaneSurface(cintfa_ar[i]) for i in range(0,4)]
114  sintfb_ar=[PlaneSurface(cintfb_ar[i]) for i in range(0,4)]  sintfb_ar=[PlaneSurface(cintfb_ar[i]) for i in range(0,4)]
115    
116  vintfa=Volume(SurfaceLoop(stop,*tuple(sintfa_ar+sintf_ar)))  sintf=PlaneSurface(CurveLoop(*tuple(linhe_ar)))
117  vintfb=Volume(SurfaceLoop(sbot,*tuple(sintfb_ar+sintf_ar)))  
118    vintfa=Volume(SurfaceLoop(stop,-sintf,*tuple(sintfa_ar)))
119    vintfb=Volume(SurfaceLoop(sbot,sintf,*tuple(sintfb_ar)))
120    
121  # Create the volume.  # Create the volume.
122  #sloop=SurfaceLoop(stop,sbot,*tuple(sintfa_ar+sintfb_ar))  #sloop=SurfaceLoop(stop,sbot,*tuple(sintfa_ar+sintfb_ar))
# Line 206  vintfb=Volume(SurfaceLoop(sbot,*tuple(si Line 125  vintfb=Volume(SurfaceLoop(sbot,*tuple(si
125    
126  #############################################EXPORTING MESH FOR ESCRIPT  #############################################EXPORTING MESH FOR ESCRIPT
127  # Create a Design which can make the mesh  # Create a Design which can make the mesh
128  d=Design(dim=3, element_size=200.0*m)  d=Design(dim=3, element_size=5.0, order=2)
 # Add the subdomains and flux boundaries.  
 #d.addItems(stop,sbot)  
 #d.addItems(PropertySet('intf',sintf))  
   
 #d.addItems(PropertySet('vintfa',vintfa))  
 #d.addItems(PropertySet('vintfb',vintfb))  
 #d.addItems(PropertySet('l1',linte_ar[0]))  
 d.addItems(*tuple(sintfb_ar))  
 d.addItems(sbot)  
 d.addItems(*tuple(sintf_ar))  
129    
130  d.setScriptFileName(os.path.join(save_path,"example09m.geo"))  d.addItems(PropertySet('vintfa',vintfa))
131    d.addItems(PropertySet('vintfb',vintfb))
132    d.addItems(PropertySet('stop',stop))
133    d.addItems(PropertySet('sbpt',sbot))
134    
135  #d.setMeshFileName(os.path.join(save_path,"example09m.msh"))  d.setScriptFileName(os.path.join(save_path,"example09m.geo"))
136    d.setMeshFileName(os.path.join(save_path,"example09m.msh"))
137  #  #
138  #  make the finley domain:  #  make the finley domain:
139  #  #
140  domain=MakeDomain(d)  domain=MakeDomain(d)
141  # Create a file that can be read back in to python with  # Create a file that can be read back in to python with
142  # mesh=ReadMesh(fileName)  # mesh=ReadMesh(fileName)
143  #domain.write(os.path.join(save_path,"example09m.fly"))  domain.write(os.path.join(save_path,"example09m.fly"))
144    
145    from esys.pycad import layer_cake
146    intfaces=[10,30,50,55,80,100,200,250,400,500]
147    
148    # Specify the domain.
149    domaindes=Design(dim=3,element_size=element_size,order=2)
150    cmplx_domain=layer_cake.LayerCake(domaindes,xwidth,ywidth,intfaces)
151    cmplx_domain.setScriptFileName(os.path.join(save_path,"example09lc.geo"))
152    cmplx_domain.setMeshFileName(os.path.join(save_path,"example09lc.msh"))
153    dcmplx=MakeDomain(cmplx_domain)
154    dcmplx.write(os.path.join(save_path,"example09lc.fly"))
155    
156    
157    
158    

Legend:
Removed from v.3068  
changed lines
  Added in v.3069

  ViewVC Help
Powered by ViewVC 1.1.26