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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 3774 - (show annotations)
Wed Jan 18 06:29:34 2012 UTC (7 years, 8 months ago) by jfenwick
File MIME type: text/x-python
File size: 5670 byte(s)
dudley, pasowrap, pycad

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.dudley 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.dudley 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.getDomain().setTagMap('tag1',1)
157 fs.setTags('tag1',m)
158 fs.setTags(100,1-m)
159 dom.write(os.path.join(MESH_DIRECTORY,filename))
160 # saveVTK(os.path.join(MESH_DIRECTORY,filename+".vtu"),dom)

  ViewVC Help
Powered by ViewVC 1.1.26