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

Revision 2914 - (show annotations)
Wed Feb 3 05:38:12 2010 UTC (9 years, 6 months ago) by gross
File MIME type: text/x-python
File size: 4398 byte(s)
missing import added.
 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 from esys.escript import mkDir 39 40 #used to construct polygons for plotting from pycad 41 from cblib1 import getLoopCoords 42 43 44 45 #################################################ESTABLISHING VARIABLES 46 #set modal to 1 for a syncline or -1 for an anticline structural 47 #configuration 48 modal=-1 49 50 # the folder to put our outputs in, leave blank "" for script path - 51 # note this folder path must exist to work 52 save_path= os.path.join("data","heatrefrac001") 53 mkDir(save_path) 54 55 ################################################ESTABLISHING PARAMETERS 56 #Model Parameters 57 width=5000.0 #width of model 58 depth=-6000.0 #depth of model 59 sspl=51 #number of discrete points in spline 60 dsp=width/(sspl-1) #dx of spline steps for width 61 dep_sp=2500.0 #avg depth of spline 62 h_sp=1500.0 #heigh of spline 63 orit=-1.0 #orientation of spline 1.0=>up -1.0=>down 64 65 ####################################################DOMAIN CONSTRUCTION 66 # Domain Corners 67 p0=Point(0.0, 0.0, 0.0) 68 p1=Point(0.0, depth, 0.0) 69 p2=Point(width, depth, 0.0) 70 p3=Point(width, 0.0, 0.0) 71 # Join corners in anti-clockwise manner. 72 l01=Line(p0, p1) 73 l12=Line(p1, p2) 74 l23=Line(p2, p3) 75 l30=Line(p3, p0) 76 # Join line segments to create domain boundary. 77 c=CurveLoop(l01, l12, l23, l30) 78 # Generate Material Boundary 79 x=[ Point(i*dsp\ 80 ,-dep_sp+modal*orit*h_sp*cos(pi*i*dsp/dep_sp+pi))\ 81 for i in range(0,sspl)\ 82 ] 83 mysp = Spline(*tuple(x)) 84 # Start and end of material boundary. 85 x1=Spline.getStartPoint(mysp) 86 x2=Spline.getEndPoint(mysp) 87 88 # Create TOP BLOCK 89 # lines 90 tbl1=Line(p0,x1) 91 tbl2=mysp 92 tbl3=Line(x2,p3) 93 # curve 94 tblockloop = CurveLoop(tbl1,tbl2,tbl3,l30) 95 # surface 96 tblock = PlaneSurface(tblockloop) 97 # polygon for python plotting 98 tpg = getLoopCoords(tblockloop) 99 np.savetxt(os.path.join(save_path,"toppg"),tpg,delimiter=" ") 100 101 # Create BOTTOM BLOCK 102 # lines 103 bbl1=Line(x1,p1) 104 bbl3=Line(p2,x2) 105 bbl4=-mysp 106 # curve 107 bblockloop = CurveLoop(bbl1,l12,bbl3,bbl4) 108 # surface 109 bblock = PlaneSurface(bblockloop) 110 111 #clockwise check as splines must be set as polygons in the point order 112 #they were created. Otherwise get a line across plot. 113 bblockloop2=CurveLoop(mysp,Line(x2,p2),Line(p2,p1),Line(p1,x1)) 114 bpg = getLoopCoords(bblockloop2) 115 np.savetxt(os.path.join(save_path,"botpg"),bpg,delimiter=" ") 116 117 #############################################EXPORTING MESH FOR ESCRIPT 118 # Create a Design which can make the mesh 119 d=Design(dim=2, element_size=200) 120 # Add the subdomains and flux boundaries. 121 d.addItems(PropertySet("top",tblock),PropertySet("bottom",bblock),\ 122 PropertySet("linebottom",l12)) 123 # Create the geometry, mesh and Escript domain 124 d.setScriptFileName(os.path.join(save_path,\ 125 "heatrefraction_mesh001.geo")) 126 127 d.setMeshFileName(os.path.join(save_path,"heatrefraction_mesh001.msh")) 128 domain=MakeDomain(d, integrationOrder=-1, reducedIntegrationOrder=-1,\ 129 optimizeLabeling=True) 130 # Create a file that can be read back in to python with 131 # mesh=ReadMesh(fileName) 132 domain.write(os.path.join(save_path,"heatrefraction_mesh001.fly")) 133 134

 ViewVC Help Powered by ViewVC 1.1.26