/[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 3346 - (show annotations)
Fri Nov 12 01:19:02 2010 UTC (9 years, 2 months ago) by caltinay
File MIME type: text/x-python
File size: 5698 byte(s)
Replaced usage of esys.escript.util.saveVTK by weipa.saveVTK in all python
scripts.

1 ########################################################
2 #
3 # Copyright (c) 2003-2010 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-2010 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 from esys.weipa import saveVTK
49 import os
50
51 def getMesh(NE_X, NE_Y, t,d,o,fullOrder,r,l_X):
52 if t == "Hex":
53 if d == 2:
54 dom=Rectangle(n0=NE_X, n1=NE_Y, l0=l_X,order=o, useFullElementOrder=fullOrder,useElementsOnFace=r,optimize=True)
55 else:
56 Brick()
57 dom=Brick(n0=NE_X, n1=NE_Y, n2=NE_Y,l0=l_X,order=o, useFullElementOrder=fullOrder,useElementsOnFace=r,optimize=True)
58 else:
59 des=Design(dim=d, order=o, element_size =min(l_X/max(3,NE_X),1./max(3,NE_Y)), keep_files=True)
60 des.setScriptFileName("tester.geo")
61 if d == 2:
62 p0=Point(0.,0.)
63 p1=Point(l_X,0.)
64 p2=Point(l_X,1.)
65 p3=Point(0.,1.)
66 l01=Line(p0,p1)
67 l12=Line(p1,p2)
68 l23=Line(p2,p3)
69 l30=Line(p3,p0)
70 s=PlaneSurface(CurveLoop(l01,l12,l23,l30))
71 des.addItems(s,l01,l12,l23,l30)
72 else:
73 p000=Point( 0.,0.,0.)
74 p100=Point(l_X,0.,0.)
75 p010=Point(0.,1.,0.)
76 p110=Point(l_X,1.,0.)
77 p001=Point(0.,0.,1.)
78 p101=Point(l_X,0.,1.)
79 p011=Point(0.,1.,1.)
80 p111=Point(l_X,1.,1.)
81
82 l10=Line(p000,p100)
83 l20=Line(p100,p110)
84 l30=Line(p110,p010)
85 l40=Line(p010,p000)
86
87 l11=Line(p000,p001)
88 l21=Line(p100,p101)
89 l31=Line(p110,p111)
90 l41=Line(p010,p011)
91
92 l12=Line(p001,p101)
93 l22=Line(p101,p111)
94 l32=Line(p111,p011)
95 l42=Line(p011,p001)
96
97 bottom=PlaneSurface(-CurveLoop(l10,l20,l30,l40))
98 top=PlaneSurface(CurveLoop(l12,l22,l32,l42))
99
100 front=PlaneSurface(CurveLoop(l10,l21,-l12,-l11))
101 back=PlaneSurface(CurveLoop(l30,l41,-l32,-l31))
102
103 left=PlaneSurface(CurveLoop(l11,-l42,-l41,l40))
104 right=PlaneSurface(CurveLoop(-l21,l20,l31,-l22))
105
106 vol=Volume(SurfaceLoop(bottom,top,front,back,left,right))
107 des.addItems(vol)
108
109 dom=MakeDomain(des)
110 return dom
111
112 # test if out put dir exists:
113 if not os.path.isdir(MESH_DIRECTORY): os.mkdir(MESH_DIRECTORY)
114
115 for d in [2,3]:
116 for o in ["1", "2", "2F"]:
117 for t in [ "Hex", "Tet" ]:
118 for c in ["Yes", "No"]:
119 for r in ["Yes", "No"]:
120 filename="mesh_D%s_o%s_T%s_Contacts%s_Rich%s.fly"%(d,o,t,c,r)
121 # certain cases are excluded:
122 if ( o == "2F" and t == "Tet" ) or \
123 ( t == "Tet" and r == "Yes" ) or \
124 ( o == "2F" and r == "Yes" ) :
125 pass
126 else:
127 print "generate file ",filename
128 if c == "Yes":
129 NE_X=int(NE**(1./d)/2+0.5)
130 NE_Y=int(NE**(1./d)+0.5)
131 else:
132 NE_X=int(NE**(1./d)+0.5)
133 NE_Y=NE_X
134 print filename
135 print "generating ",NE_X, NE_Y
136 if o == "2":
137 o2=2
138 full=False
139 elif o == "2F":
140 o2=2
141 full=True
142 else:
143 o2=1
144 full=False
145 if c == "Yes":
146 dom1=getMesh(NE_X, NE_Y,t,d,o2,full,r=="Yes",0.5)
147 dom2=getMesh(NE_X, NE_Y,t,d,o2,full,r=="Yes",0.5)
148 x=dom2.getX().copy()
149 x[0]=1.-x[0]
150 dom2.setX(x)
151 dom=JoinFaces([dom1,dom2])
152 else:
153 dom=getMesh(NE_X, NE_Y,t,d,o2,full,r=="Yes",1.)
154 # set tags:
155 for fs in [ContinuousFunction(dom), Function(dom), FunctionOnBoundary(dom), FunctionOnContactOne(dom)]:
156 m=whereNegative(fs.getX()[d-1]-CUT)
157 fs.getDomain().setTagMap('tag1',1)
158 fs.setTags('tag1',m)
159 fs.setTags(100,1-m)
160 dom.write(os.path.join(MESH_DIRECTORY,filename))
161 # saveVTK(os.path.join(MESH_DIRECTORY,filename+".vtu"),dom)

  ViewVC Help
Powered by ViewVC 1.1.26