/[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 2841 - (hide annotations)
Thu Jan 14 01:02:56 2010 UTC (9 years, 8 months ago) by jfenwick
File MIME type: text/x-python
File size: 5667 byte(s)
Unit tests for 2838 and adding setTags variant which takes a string instead of a tag num.


1 jfenwick 1851 ########################################################
2     #
3 jfenwick 2548 # Copyright (c) 2003-2009 by University of Queensland
4 jfenwick 1851 # 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 jfenwick 2549 __copyright__="""Copyright (c) 2003-2009 by University of Queensland
14 jfenwick 1851 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 jfenwick 2344 __url__="https://launchpad.net/escript-finley"
20 jfenwick 1851 """
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 jfenwick 2841 fs.getDomain().setTagMap('tag1',1)
157     fs.setTags('tag1',m)
158 jfenwick 1851 fs.setTags(100,1-m)
159     dom.write(os.path.join(MESH_DIRECTORY,filename))
160 caltinay 2534 # saveVTK(os.path.join(MESH_DIRECTORY,filename+".vtu"),dom)

  ViewVC Help
Powered by ViewVC 1.1.26