← Back to branch summary

~reducedmodelling/fluidity/ROM_Non-intrusive-ann

« back to all changes in this revision

Viewing changes to tests/turbine_flux_dg_2d/mesh/scripts/triangle.py

  • Committer: sf1409
  • Date: 2010-05-10 15:46:37 UTC
  • Revision ID: svn-v4:5bf5533e-7014-46e3-b1bb-cce4b9d03719:trunk:1140
First of four flux turbine unleashed:

This is a 2d flux turbine testcase.
It uses dg flux terms to simulate the turbine. This is a very basic testcase where the turbine is basically a opened sluice gate. The result is compared with another simulation where no turbine is applied. Ideally, the results should be the same. Due to the different meshes, the detectors will have slightly different values.


Show diffs side-by-side

added added

removed removed

Lines of Context:
tests/turbine_flux_dg_2d/mesh/scripts/triangle.py
 
1
#!/usr/bin/python
 
2
import sys
 
3
from sets import Set
 
4
 
 
5
edges=[]
 
6
surface_ids=[] # surface ids of the edges
 
7
nodes=[]
 
8
eles=[]
 
9
region_ids=[] # region ids of the elements
 
10
edge_attribute1=[] # first edge attribute 
 
11
dim=-1
 
12
debug=0
 
13
 
 
14
 
 
15
def read_nodefile(filename):
 
16
   f = open(filename, 'r')
 
17
   global nodes, dim
 
18
   [nodes, dim]=pnode(f)
 
19
   f.close()
 
20
 
 
21
def read_elefile(filename):
 
22
   f = open(filename, 'r')
 
23
   global eles
 
24
   eles=peles(f)
 
25
   f.close()
 
26
 
 
27
def read_edgefile(filename):
 
28
   f = open(filename, 'r')
 
29
   global edges, surface_ids, edge_attribute1
 
30
   [edges, surface_ids, edge_attribute1]=pedges(f)
 
31
   f.close()
 
32
 
 
33
def save_nodefile(nodes, dim, filename):
 
34
   f = open(filename, 'w')
 
35
 
 
36
   f.write(`len(nodes)`+' '+`dim`+' 0 0\n')
 
37
   k=1
 
38
   for node in nodes:
 
39
        f.write(`k`+' ' + ''.join([`n`+' ' for n in node])+'\n')
 
40
        k=k+1
 
41
   f.write('# Created with triangle.py')
 
42
   f.close()
 
43
 
 
44
def save_elefile(eles, par1, filename):
 
45
   f = open(filename, 'w')
 
46
   f.write(`len(eles)`+' '+`len(eles[0])`+' 1\n')
 
47
   if len(par1)!=len(eles):
 
48
        print "Error in save_elefile(): Length of element list is ", len(eles), "  but length of region_id list is: ", len(par1)
 
49
        exit()
 
50
   k=1
 
51
   for ele in eles:
 
52
        f.write(`k`+' ' + ''.join([`e`+' ' for e in ele])+`par1[k-1]`+'\n')
 
53
        k=k+1
 
54
   f.write('# Created with triangle.py')
 
55
   f.close()
 
56
 
 
57
 
 
58
def save_edgefile(edges, atr1, filename):
 
59
   f = open(filename, 'w')
 
60
   f.write(`len(edges)`+' 1\n')
 
61
   if len(atr1)!=len(edges):
 
62
        print "Error in save_edgefile(): Length of edge list is ", len(edges), "  but length of atr1 list is: ", len(par1)
 
63
        exit()
 
64
   k=1
 
65
   for edge in edges:
 
66
        f.write(`k`+' ' + ''.join([`e`+' ' for e in edge])+`atr1[k-1]`+'\n')
 
67
        k=k+1
 
68
   f.write('# Create with triangly.py')
 
69
   f.close()
 
70
 
 
71
def save_edgefile2(edges, atr1, atr2, filename):
 
72
   f = open(filename, 'w')
 
73
   f.write(`len(edges)`+' 2\n')
 
74
   if len(atr1)!=len(edges) or len(atr2)!=len(edges):
 
75
        print "Error in save_edgefile(): Length of edge list is ", len(edges), "  but length of atr1, atr2 list is: ", len(atr1), ', ', len(atr2)
 
76
        exit()
 
77
   k=1
 
78
   for edge in edges:
 
79
        f.write(`k`+' ' + ''.join([`e`+' ' for e in edge])+`atr1[k-1]`+' '+ `atr2[k-1]`+'\n')
 
80
        k=k+1
 
81
   f.write('# Create with triangle.py')
 
82
   f.close()
 
83
 
 
84
def pnode(fnode):
 
85
    """ Parse the .node file
 
86
    and returns the correspondind
 
87
    list of the nodes
 
88
    """
 
89
 
 
90
    nodes = []
 
91
    count = 0
 
92
    dim = 0
 
93
    nbnodes = 0
 
94
 
 
95
    for line in fnode:
 
96
 
 
97
        if line[0]=='#':
 
98
            continue
 
99
        else:
 
100
            parts = line.split()
 
101
            if len(parts)==0:
 
102
                continue
 
103
            if count == 0:
 
104
                nbnodes = int(parts[0])
 
105
                dim = int(parts[1])
 
106
                count += 1
 
107
                if dim not in [2,3]:
 
108
                    print "Error dim in .node file!"
 
109
 
 
110
            else:
 
111
                if dim == 3:
 
112
                    nodes.append(([float(x) for x in parts[1:]]))
 
113
                else:
 
114
                    nodes.append(([float(x) for x in parts[1:3]]))
 
115
    if (len(nodes) != nbnodes):
 
116
        print "Error: actual number of nodes does not corresponds " + str(len(nodes)), nbnodes
 
117
        sys.exit(2)
 
118
                                 
 
119
    return nodes, dim
 
120
 
 
121
def pedges(fedge):
 
122
    """ Parse the .edge file and returns the corresponding
 
123
    list of edges
 
124
    """
 
125
 
 
126
    edges = []
 
127
    surface_ids=[]
 
128
    attribute1=[]
 
129
    count = 0
 
130
    nbedges = 0
 
131
    boundary = 0
 
132
    for line in fedge:
 
133
        if line[0]=='#':
 
134
            continue
 
135
        else:
 
136
            parts = line.split()
 
137
            if len(parts)==0:
 
138
                continue
 
139
        
 
140
            if len(parts)==0:
 
141
                continue
 
142
            if count == 0:
 
143
                nbedges = int(parts[0])
 
144
                boundary = int(parts[1])
 
145
                count += 1
 
146
                if boundary>2:
 
147
                    print "Error: too many attributes in .edge file!"
 
148
 
 
149
            else:
 
150
                 edges.append(([int(x) for x in parts[1:3]]))
 
151
                 if boundary >= 1:
 
152
                   surface_ids.append(int(parts[3]))
 
153
                 if boundary >= 2:
 
154
                   attribute1.append(int(parts[4]))
 
155
 
 
156
    if (len(edges) != nbedges):
 
157
        print "Error: actual number of edges does not corresponds " + str(len(edges)), str(nbedges)
 
158
        sys.exit(2)
 
159
 
 
160
    return [edges, surface_ids, attribute1]
 
161
 
 
162
def peles(fele):
 
163
    """ Parse the .ele file and returns the corresponding
 
164
    list of elements
 
165
    """
 
166
 
 
167
    eles = []
 
168
    count = 0
 
169
    nbeles = 0
 
170
    for line in fele:
 
171
 
 
172
        if line[0]=='#':
 
173
            continue
 
174
        else:
 
175
            parts = line.split()
 
176
 
 
177
            if count == 0:
 
178
                nbeles = int(parts[0])
 
179
                nodesperele = int(parts[1])
 
180
                region = int(parts[2])
 
181
                count += 1
 
182
 
 
183
            else:
 
184
                eles.append(([int(x) for x in parts[1:nodesperele+1]]))
 
185
                if region==1:
 
186
                        region_ids.append(int(parts[len(parts)-1]))
 
187
 
 
188
    if (len(eles) != nbeles):
 
189
        print "Error: actual number of eles does not corresponds " + str(len(eles)), str(nbeles)
 
190
        sys.exit(2)
 
191
 
 
192
    return eles
 
193
 
 
194
# Returns array of edge indices whose boundary marker ids equal bid
 
195
def edges_with_bid(bid):
 
196
   global edges
 
197
   ret=[]
 
198
   
 
199
   for i in range(0,len(surface_ids)):
 
200
        if surface_ids[i] == bid:
 
201
                ret.append(i)
 
202
   return ret
 
203
 
 
204
# Returns array of ele indices which have the given edge id
 
205
def ele_with_edgeids(edgeids):
 
206
   global edges, eles, nodes
 
207
   ret=[]
 
208
   nodeids=edges[edgeids][0:len(edges[edgeids])] # get nodeids of edgeid
 
209
   # loop through all elements and check if they include all these nodes
 
210
   i=0
 
211
   for ele in eles:
 
212
        if Set(nodeids).issubset(Set(ele[0:len(ele)])):
 
213
                ret.append(i)
 
214
        i=i+1
 
215
   if len(ret)>2 :
 
216
        print "There are ", len(ret), " elements adjacent to face with id ", edgeids, ". Only 2 were expected."
 
217
        exit()
 
218
   return ret
 
219
 
 
220
# Returns array of ele indices which have the given node ids
 
221
def ele_with_nodeids(nodeids):
 
222
   global eles, nodes
 
223
   ret=[]
 
224
   # loop through all elements and check if they include all these nodes
 
225
   i=0
 
226
   for ele in eles:
 
227
        if Set(nodeids).issubset(Set(ele[0:len(ele)])):
 
228
                ret.append(i)
 
229
        i=i+1
 
230
   if len(ret)!=2 :
 
231
        print "There are ", len(ret), " elements adjacent to nodes with ids ", nodeids, ". Only 2 were expected."
 
232
        exit()
 
233
   return ret
 
234
 
 
235
# returns the edges with nodeids
 
236
def edgeid_from_nodeids(nodeids):
 
237
        global edges
 
238
        ret=[]
 
239
        i=0
 
240
        for edge in edges:
 
241
                if Set(nodeids).issubset(Set(edge[0:len(edge)])):
 
242
                        ret.append(i)
 
243
                i=i+1
 
244
        return ret
 
245
 
 
246
# returns the edge with boundary_id and wit node nodeid but not with id edgeid
 
247
def get_partner_edge(edgein, nodeid, boundary_id):
 
248
        global edges
 
249
        cands=edges_with_bid(boundary_id)
 
250
        for edge in cands:
 
251
                if edge!=edgein and (nodeid in edges[edge]):
 
252
                        return edge
 
253
 
 
254
def get_boundaryid_from_edgeid(edgeid):
 
255
        global surface_ids
 
256
        return surface_ids[edgeid]
 
257
 
 
258
def has_ele_edge_on_boundaryid(elein, nodeid, boundary_id):
 
259
        global eles, debug
 
260
        for n in eles[elein]:
 
261
                if n==nodeid:
 
262
                        continue
 
263
                edge=edgeid_from_nodeids([n,nodeid])
 
264
                if edge!=[] and get_boundaryid_from_edgeid(edge[0])==boundary_id:
 
265
                        if debug>3:
 
266
                                print "Found boundary edge with id", edgeid_from_nodeids([n,nodeid])[0]+1
 
267
                        return True
 
268
        return False
 
269
                
 
270
def get_eles_on_ele_side(elein, nodeid, edgein, boundary_id):
 
271
        global eles, edges
 
272
        if not has_ele_edge_on_boundaryid(elein, nodeid, boundary_id):
 
273
                print "Error in get_eles_on_ele_side: Given Element has no edge with secified boundary id"
 
274
                sys.exit()
 
275
        ret=[elein]
 
276
        thirdnode=Set(eles[elein]).difference(Set(edges[edgein])).pop()
 
277
        while True: 
 
278
                if debug>3:
 
279
                        print "Found third node of current element wit index: ", thirdnode
 
280
                        print "Looking for neighbour with nodes ", thirdnode, nodeid
 
281
                partner_eleins=ele_with_nodeids([thirdnode, nodeid])
 
282
                if partner_eleins[0]==elein:
 
283
                        partner_elein=partner_eleins[1]
 
284
                else:
 
285
                        partner_elein=partner_eleins[0]
 
286
                if debug>3:
 
287
                        print "Found neighbour with eleid", partner_elein+1
 
288
                ret.append(partner_elein)
 
289
                if has_ele_edge_on_boundaryid(partner_elein, nodeid, boundary_id):
 
290
                        if debug>3:
 
291
                                print "Reached the end of the element loop"
 
292
                        return ret
 
293
                if debug>3:
 
294
                        print "New elemnt has no desired boundary edge, loop again"
 
295
                elein=partner_elein
 
296
                thirdnode=Set(eles[elein]).difference(Set([nodeid, thirdnode])).pop()
 
297
 
 
298
 
 
299
# deletes the node with id nodeid and updates the elelist and edgelist 
 
300
def delete_nodeid(nodeid):
 
301
        global eles, edges, nodes
 
302
        for index, ele in enumerate(eles):
 
303
                if nodeid in ele:
 
304
                        print "Error in delete_nodeid(", nodeid, "). This node is still in use in element ", index+1, ": ", ele, "."
 
305
                        exit()
 
306
                eles[index]=[x if x<nodeid else x-1 for x in eles[index]]
 
307
        for index, edge in enumerate(edges):
 
308
                if nodeid in edge:
 
309
                        print "Error in delete_nodeid(", nodeid, "). This node is still in use in edge ", index+1, "."
 
310
                        exit()
 
311
                edges[index]=[x if x<nodeid else x-1 for x in edges[index]]
 
312
        nodes.pop(nodeid-1)
 
313
 
 
314
# deletes the ele with id eleid 
 
315
def delete_eleid(eleid):
 
316
        global eles, region_ids
 
317
        eles.pop(eleid-1)
 
318
        region_ids.pop(eleid-1)
 
319
 
 
320
# deletes the edge with id edgeid
 
321
def delete_edgeid(edgeid):
 
322
        global edges, surface_ids
 
323
        edges.pop(edgeid-1)
 
324
        surface_ids.pop(nodeid-1)
 
325
        edge_attribute1.pop(nodeid-1)