# Contents of /trunk/doc/examples/cookbook/heatrefraction_mesher001.py

Revision 2904 - (show annotations)
Tue Feb 2 04:21:52 2010 UTC (9 years, 2 months ago) by gross
File MIME type: text/x-python
File size: 4367 byte(s)
```more revision on cookbook: example 1 and 2 have been swapped.
```
 1 2 ######################################################## 3 # 4 # Copyright (c) 2009-2010 by University of Queensland 5 # Earth Systems Science Computational Center (ESSCC) 6 7 # 8 # Primary Business: Queensland, Australia 9 # Licensed under the Open Software License version 3.0 10 11 # 12 ######################################################## 13 14 __copyright__="""Copyright (c) 2009-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 20 __url__= 21 22 """ 23 Author: Antony Hallam antony.hallam@uqconnect.edu.au 24 """ 25 26 ############################################################FILE HEADER 27 # heatrefraction_mesher001.py 28 # Create either a 2D syncline or anticline model using pycad meshing 29 # tools. 30 31 #######################################################EXTERNAL MODULES 32 from esys.pycad import * #domain constructor 33 from esys.pycad.gmsh import Design #Finite Element meshing package 34 from esys.finley import MakeDomain #Converter for escript 35 import os #file path tool 36 import numpy as np #numerial python for arrays 37 from math import * # math package 38 39 #used to construct polygons for plotting from pycad 40 from cblib1 import getLoopCoords 41 42 43 44 #################################################ESTABLISHING VARIABLES 45 #set modal to 1 for a syncline or -1 for an anticline structural 46 #configuration 47 modal=-1 48 49 # the folder to put our outputs in, leave blank "" for script path - 50 # note this folder path must exist to work 51 save_path= os.path.join("data","heatrefrac001") 52 mkdir(save_path) 53 54 ################################################ESTABLISHING PARAMETERS 55 #Model Parameters 56 width=5000.0 #width of model 57 depth=-6000.0 #depth of model 58 sspl=51 #number of discrete points in spline 59 dsp=width/(sspl-1) #dx of spline steps for width 60 dep_sp=2500.0 #avg depth of spline 61 h_sp=1500.0 #heigh of spline 62 orit=-1.0 #orientation of spline 1.0=>up -1.0=>down 63 64 ####################################################DOMAIN CONSTRUCTION 65 # Domain Corners 66 p0=Point(0.0, 0.0, 0.0) 67 p1=Point(0.0, depth, 0.0) 68 p2=Point(width, depth, 0.0) 69 p3=Point(width, 0.0, 0.0) 70 # Join corners in anti-clockwise manner. 71 l01=Line(p0, p1) 72 l12=Line(p1, p2) 73 l23=Line(p2, p3) 74 l30=Line(p3, p0) 75 # Join line segments to create domain boundary. 76 c=CurveLoop(l01, l12, l23, l30) 77 # Generate Material Boundary 78 x=[ Point(i*dsp\ 79 ,-dep_sp+modal*orit*h_sp*cos(pi*i*dsp/dep_sp+pi))\ 80 for i in range(0,sspl)\ 81 ] 82 mysp = Spline(*tuple(x)) 83 # Start and end of material boundary. 84 x1=Spline.getStartPoint(mysp) 85 x2=Spline.getEndPoint(mysp) 86 87 # Create TOP BLOCK 88 # lines 89 tbl1=Line(p0,x1) 90 tbl2=mysp 91 tbl3=Line(x2,p3) 92 # curve 93 tblockloop = CurveLoop(tbl1,tbl2,tbl3,l30) 94 # surface 95 tblock = PlaneSurface(tblockloop) 96 # polygon for python plotting 97 tpg = getLoopCoords(tblockloop) 98 np.savetxt(os.path.join(save_path,"toppg"),tpg,delimiter=" ") 99 100 # Create BOTTOM BLOCK 101 # lines 102 bbl1=Line(x1,p1) 103 bbl3=Line(p2,x2) 104 bbl4=-mysp 105 # curve 106 bblockloop = CurveLoop(bbl1,l12,bbl3,bbl4) 107 # surface 108 bblock = PlaneSurface(bblockloop) 109 110 #clockwise check as splines must be set as polygons in the point order 111 #they were created. Otherwise get a line across plot. 112 bblockloop2=CurveLoop(mysp,Line(x2,p2),Line(p2,p1),Line(p1,x1)) 113 bpg = getLoopCoords(bblockloop2) 114 np.savetxt(os.path.join(save_path,"botpg"),bpg,delimiter=" ") 115 116 #############################################EXPORTING MESH FOR ESCRIPT 117 # Create a Design which can make the mesh 118 d=Design(dim=2, element_size=200) 119 # Add the subdomains and flux boundaries. 120 d.addItems(PropertySet("top",tblock),PropertySet("bottom",bblock),\ 121 PropertySet("linebottom",l12)) 122 # Create the geometry, mesh and Escript domain 123 d.setScriptFileName(os.path.join(save_path,\ 124 "heatrefraction_mesh001.geo")) 125 126 d.setMeshFileName(os.path.join(save_path,"heatrefraction_mesh001.msh")) 127 domain=MakeDomain(d, integrationOrder=-1, reducedIntegrationOrder=-1,\ 128 optimizeLabeling=True) 129 # Create a file that can be read back in to python with 130 # mesh=ReadMesh(fileName) 131 domain.write(os.path.join(save_path,"heatrefraction_mesh001.fly")) 132 133