/[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 1851 - (hide annotations)
Mon Oct 6 03:16:43 2008 UTC (11 years ago) by jfenwick
Original Path: branches/more_shared_ptrs_from_1812/finley/test/python/generate_meshes.py
File MIME type: text/x-python
File size: 5623 byte(s)
Branch commit.
Added files while the merge op did not add.
Modified shake59 config for non-ken users.

1 jfenwick 1851 ########################################################
2     #
3     # Copyright (c) 2003-2008 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-2008 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__="http://www.uq.edu.au/esscc/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     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.setTags(1,m)
157     fs.setTags(100,1-m)
158     dom.write(os.path.join(MESH_DIRECTORY,filename))
159     # saveVTK(os.path.join(MESH_DIRECTORY,filename+".xml"),dom)

  ViewVC Help
Powered by ViewVC 1.1.26