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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 6651 - (hide annotations)
Wed Feb 7 02:12:08 2018 UTC (19 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 jfenwick 3981 ##############################################################################
2 jfenwick 1851 #
3 jfenwick 6651 # Copyright (c) 2003-2018 by The University of Queensland
4 jfenwick 3981 # http://www.uq.edu.au
5 jfenwick 1851 #
6     # Primary Business: Queensland, Australia
7 jfenwick 6112 # Licensed under the Apache License, version 2.0
8     # http://www.apache.org/licenses/LICENSE-2.0
9 jfenwick 1851 #
10 jfenwick 3981 # Development until 2012 by Earth Systems Science Computational Center (ESSCC)
11 jfenwick 4657 # Development 2012-2013 by School of Earth Sciences
12     # Development from 2014 by Centre for Geoscience Computing (GeoComp)
13 jfenwick 3981 #
14     ##############################################################################
15 jfenwick 1851
16 sshaw 5706 from __future__ import print_function, division
17    
18 jfenwick 6651 __copyright__="""Copyright (c) 2003-2018 by The University of Queensland
19 jfenwick 3981 http://www.uq.edu.au
20 jfenwick 1851 Primary Business: Queensland, Australia"""
21 jfenwick 6112 __license__="""Licensed under the Apache License, version 2.0
22     http://www.apache.org/licenses/LICENSE-2.0"""
23 jfenwick 2344 __url__="https://launchpad.net/escript-finley"
24 jfenwick 1851 """
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 caltinay 3346 from esys.weipa import saveVTK
53 jfenwick 1851 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 jfenwick 3772 print("generate file ",filename)
132 jfenwick 1851 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 jfenwick 3772 print(filename)
139     print("generating ",NE_X, NE_Y)
140 jfenwick 1851 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 sshaw 4576 fs.getDomain().setTagMap('tag1',1)
162 jfenwick 2841 fs.setTags('tag1',m)
163 jfenwick 1851 fs.setTags(100,1-m)
164     dom.write(os.path.join(MESH_DIRECTORY,filename))
165 caltinay 2534 # saveVTK(os.path.join(MESH_DIRECTORY,filename+".vtu"),dom)

  ViewVC Help
Powered by ViewVC 1.1.26