/[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 2534 - (show annotations)
Thu Jul 16 06:49:19 2009 UTC (10 years, 2 months ago) by caltinay
File MIME type: text/x-python
File size: 5618 byte(s)
Changed examples, tests and tutorials to save VTK files as .vtu instead .xml.
Visit doesn't know what to do with xml's and vtu is the proper extension
anyway.

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

  ViewVC Help
Powered by ViewVC 1.1.26