← Back to branch summary

~vcs-imports/escript-finley/trunk

« back to all changes in this revision

Viewing changes to dudley/test/python/generate_meshes.py

  • Committer: jfenwick
  • Date: 2010-10-11 01:48:14 UTC
  • Revision ID: svn-v4:77569008-7704-0410-b7a0-a92fef0b09fd:trunk:3259
Merging dudley and scons updates from branches

Show diffs side-by-side

added added

removed removed

Lines of Context:
dudley/test/python/generate_meshes.py
 
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)