/[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 6651 - (show annotations)
Wed Feb 7 02:12:08 2018 UTC (20 months, 1 week ago) by jfenwick
File MIME type: text/x-python
File size: 5882 byte(s)
Make everyone sad by touching all the files

Copyright dates update

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

  ViewVC Help
Powered by ViewVC 1.1.26