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

Contents of /trunk/finley/test/python/generate_meshes.py

Parent Directory Parent Directory | Revision Log Revision Log


Revision 3981 - (show annotations)
Fri Sep 21 02:47:54 2012 UTC (7 years, 2 months ago) by jfenwick
File MIME type: text/x-python
File size: 5763 byte(s)
First pass of updating copyright notices
1 ##############################################################################
2 #
3 # Copyright (c) 2003-2012 by University of Queensland
4 # http://www.uq.edu.au
5 #
6 # Primary Business: Queensland, Australia
7 # Licensed under the Open Software License version 3.0
8 # http://www.opensource.org/licenses/osl-3.0.php
9 #
10 # Development until 2012 by Earth Systems Science Computational Center (ESSCC)
11 # Development since 2012 by School of Earth Sciences
12 #
13 ##############################################################################
14
15 __copyright__="""Copyright (c) 2003-2012 by University of Queensland
16 http://www.uq.edu.au
17 Primary Business: Queensland, Australia"""
18 __license__="""Licensed under the Open Software License version 3.0
19 http://www.opensource.org/licenses/osl-3.0.php"""
20 __url__="https://launchpad.net/escript-finley"
21 """
22 this script generates a variety of meshes in the fly format on a unit square with a
23 fixed number of elements NE (approximatively).
24
25 This script will not run under MPI.
26
27 generated file names are:
28
29 mesh_D<d>_o<O>_T<t>_Contacts<c>_Rich<r>.fly
30
31 wehere <d>= spatial dimension 2,3
32 <o>= order 1,2, 2F
33 <t>= element type = hex, tet
34 <c>= Yes if contact elements are present, No otherwise.
35 <r>= YEs is rich surface elements are used, No otherwise
36
37 any item with x[d]<CUT will be tagged with 1 the rest will be marked 100.
38 """
39
40 MESH_DIRECTORY="./tmp_meshes"
41 NE=10
42 CUT=0.5
43
44 from esys.escript import *
45 from esys.finley import Rectangle, Brick, Merge, JoinFaces
46 from esys.pycad import Point, Line,PlaneSurface, CurveLoop, Volume,SurfaceLoop
47 from esys.pycad.gmsh import Design
48 from esys.finley import MakeDomain
49 from esys.weipa import saveVTK
50 import os
51
52 def getMesh(NE_X, NE_Y, t,d,o,fullOrder,r,l_X):
53 if t == "Hex":
54 if d == 2:
55 dom=Rectangle(n0=NE_X, n1=NE_Y, l0=l_X,order=o, useFullElementOrder=fullOrder,useElementsOnFace=r,optimize=True)
56 else:
57 Brick()
58 dom=Brick(n0=NE_X, n1=NE_Y, n2=NE_Y,l0=l_X,order=o, useFullElementOrder=fullOrder,useElementsOnFace=r,optimize=True)
59 else:
60 des=Design(dim=d, order=o, element_size =min(l_X/max(3,NE_X),1./max(3,NE_Y)), keep_files=True)
61 des.setScriptFileName("tester.geo")
62 if d == 2:
63 p0=Point(0.,0.)
64 p1=Point(l_X,0.)
65 p2=Point(l_X,1.)
66 p3=Point(0.,1.)
67 l01=Line(p0,p1)
68 l12=Line(p1,p2)
69 l23=Line(p2,p3)
70 l30=Line(p3,p0)
71 s=PlaneSurface(CurveLoop(l01,l12,l23,l30))
72 des.addItems(s,l01,l12,l23,l30)
73 else:
74 p000=Point( 0.,0.,0.)
75 p100=Point(l_X,0.,0.)
76 p010=Point(0.,1.,0.)
77 p110=Point(l_X,1.,0.)
78 p001=Point(0.,0.,1.)
79 p101=Point(l_X,0.,1.)
80 p011=Point(0.,1.,1.)
81 p111=Point(l_X,1.,1.)
82
83 l10=Line(p000,p100)
84 l20=Line(p100,p110)
85 l30=Line(p110,p010)
86 l40=Line(p010,p000)
87
88 l11=Line(p000,p001)
89 l21=Line(p100,p101)
90 l31=Line(p110,p111)
91 l41=Line(p010,p011)
92
93 l12=Line(p001,p101)
94 l22=Line(p101,p111)
95 l32=Line(p111,p011)
96 l42=Line(p011,p001)
97
98 bottom=PlaneSurface(-CurveLoop(l10,l20,l30,l40))
99 top=PlaneSurface(CurveLoop(l12,l22,l32,l42))
100
101 front=PlaneSurface(CurveLoop(l10,l21,-l12,-l11))
102 back=PlaneSurface(CurveLoop(l30,l41,-l32,-l31))
103
104 left=PlaneSurface(CurveLoop(l11,-l42,-l41,l40))
105 right=PlaneSurface(CurveLoop(-l21,l20,l31,-l22))
106
107 vol=Volume(SurfaceLoop(bottom,top,front,back,left,right))
108 des.addItems(vol)
109
110 dom=MakeDomain(des)
111 return dom
112
113 # test if out put dir exists:
114 if not os.path.isdir(MESH_DIRECTORY): os.mkdir(MESH_DIRECTORY)
115
116 for d in [2,3]:
117 for o in ["1", "2", "2F"]:
118 for t in [ "Hex", "Tet" ]:
119 for c in ["Yes", "No"]:
120 for r in ["Yes", "No"]:
121 filename="mesh_D%s_o%s_T%s_Contacts%s_Rich%s.fly"%(d,o,t,c,r)
122 # certain cases are excluded:
123 if ( o == "2F" and t == "Tet" ) or \
124 ( t == "Tet" and r == "Yes" ) or \
125 ( o == "2F" and r == "Yes" ) :
126 pass
127 else:
128 print("generate file ",filename)
129 if c == "Yes":
130 NE_X=int(NE**(1./d)/2+0.5)
131 NE_Y=int(NE**(1./d)+0.5)
132 else:
133 NE_X=int(NE**(1./d)+0.5)
134 NE_Y=NE_X
135 print(filename)
136 print("generating ",NE_X, NE_Y)
137 if o == "2":
138 o2=2
139 full=False
140 elif o == "2F":
141 o2=2
142 full=True
143 else:
144 o2=1
145 full=False
146 if c == "Yes":
147 dom1=getMesh(NE_X, NE_Y,t,d,o2,full,r=="Yes",0.5)
148 dom2=getMesh(NE_X, NE_Y,t,d,o2,full,r=="Yes",0.5)
149 x=dom2.getX().copy()
150 x[0]=1.-x[0]
151 dom2.setX(x)
152 dom=JoinFaces([dom1,dom2])
153 else:
154 dom=getMesh(NE_X, NE_Y,t,d,o2,full,r=="Yes",1.)
155 # set tags:
156 for fs in [ContinuousFunction(dom), Function(dom), FunctionOnBoundary(dom), FunctionOnContactOne(dom)]:
157 m=whereNegative(fs.getX()[d-1]-CUT)
158 fs.getDomain().setTagMap('tag1',1)
159 fs.setTags('tag1',m)
160 fs.setTags(100,1-m)
161 dom.write(os.path.join(MESH_DIRECTORY,filename))
162 # saveVTK(os.path.join(MESH_DIRECTORY,filename+".vtu"),dom)

  ViewVC Help
Powered by ViewVC 1.1.26