source: anuga_core/source/anuga/mesh_engine/mesh_engine.py @ 4898

Last change on this file since 4898 was 4898, checked in by duncan, 16 years ago

Outputting numeric arrays from mesh_engine. they aren't used yet though.

File size: 5.7 KB
Line 
1#!/usr/bin/env python
2
3import sys
4
5from types import ListType, TupleType
6
7import anuga.mesh_engine.mesh_engine_c_layer as triang
8#import anuga.mesh_engine.list_dic as triang
9
10from Numeric import array, Float, Int32, reshape
11
12from anuga.utilities.numerical_tools import ensure_numeric
13from anuga.utilities.anuga_exceptions import ANUGAError
14   
15def generate_mesh(points=None,
16                  segments=None,holes=None,regions=None,
17                  pointatts=None,segatts=None,
18                  mode=None, dummy_test=None):
19    """
20    pointatts can be a list of lists.
21
22    generatedtriangleattributelist is used to represent tagged regions.
23    #FIXME (DSG-DSG): add comments
24    """
25    #FIXME (DSG-DSG): Catch parameters that are lists,
26    #instead of lists of lists
27    # check shape[1] is 2 etc
28
29    if points is None:
30        points = []
31
32    if segments is None:
33        segments = []
34
35    if holes is None:
36        holes = []
37       
38    if regions is None:
39        regions = []
40
41    if dummy_test is None:
42        dummy_test  = []
43       
44    try:
45        points =  ensure_numeric(points, Float)
46    except ValueError:
47        msg = 'ERROR: Inconsistent points array.'
48        raise ANUGAError, msg
49    if points.shape[1] <>2:
50        msg = 'ERROR: Bad shape points array.'
51        raise ANUGAError, msg
52
53    #print "pointatts",pointatts
54    # This is after points is numeric
55    if pointatts is None or pointatts == []:
56        pointatts = [[] for x in range(points.shape[0])]
57    try:
58        # If Int is used, instead of Int32, it fails in Linux
59        segments = ensure_numeric(segments, Int32)
60       
61    except ValueError:
62        msg = 'ERROR: Inconsistent segments array.'
63        raise ANUGAError, msg
64   
65    # This is after segments is numeric
66    if segatts is None or segatts == []:
67        segatts = [0 for x in range(segments.shape[0])]
68    try:
69        holes = ensure_numeric(holes, Float)
70    except ValueError:
71        msg = 'ERROR: Inconsistent holess array.'
72        raise ANUGAError, msg
73
74   
75    regions = add_area_tag(regions)
76    try:
77        regions = ensure_numeric(regions, Float)
78    except  (ValueError, TypeError):
79        msg = 'ERROR: Inconsistent regions array.'
80        raise ANUGAError, msg
81       
82    if not regions.shape[0] == 0 and regions.shape[1] <= 2:
83        msg = 'ERROR: Bad shape points array.'
84        raise ANUGAError, msg
85   
86    try:
87        pointatts = ensure_numeric(pointatts, Float)
88    except (ValueError, TypeError):
89        msg = 'ERROR: Inconsistent point attributes array.'
90        raise ANUGAError, msg
91
92    if pointatts.shape[0] <> points.shape[0]:
93        msg = """ERROR: Point attributes array not the same shape as
94        point array."""
95        raise ANUGAError, msg
96    if len(pointatts.shape) == 1:
97        pointatts = reshape(pointatts,(pointatts.shape[0],1))
98   
99    try:
100        segatts = ensure_numeric(segatts, Int32)
101    except ValueError:
102        msg = 'ERROR: Inconsistent point attributes array.'
103        raise ANUGAError, msg
104    if segatts.shape[0] <> segments.shape[0]:
105        msg = """ERROR: Segment attributes array not the same shape as
106        segment array."""
107        raise ANUGAError, msg
108
109   
110    #print "mode", mode
111    if mode.find('n'):
112        #pass
113        mode = 'j' + mode
114        # j- Jettisons vertices that are not part of the final
115        #    triangulation from the output .node file (including duplicate
116        #    input vertices and vertices ``eaten'' by holes).  - output a
117        #    list of neighboring triangles
118        # EG handles lone verts!
119           
120    #print "points",points
121    #print "segments", segments
122    #print "segments.shape", segments.shape
123    #print "holes", holes
124    #print "regions", regions
125    #print "pointatts", pointatts
126    #print "segatts", segatts
127    #XSprint "mode", mode
128    #print "yeah"
129    mesh_dict, trianglelist, pointlist, pointmarkerlist, pointattributelist, triangleattributelist, segmentlist, segmentmarkerlist, neighborlist = triang.genMesh(points,segments,holes,regions,
130                          pointatts,segatts, mode, segments.flat)
131    # the values as arrays
132    mesh_dict['trianglelist'] = trianglelist
133    mesh_dict['pointlist'] = pointlist
134
135    # WARNING - pointmarkerlist IS UNTESTED
136    mesh_dict['pointmarkerlist'] = pointmarkerlist
137    mesh_dict['pointattributelist'] = pointattributelist
138    mesh_dict['triangleattributelist'] = triangleattributelist
139    mesh_dict['segmentlist'] = segmentlist
140    mesh_dict['segmentmarkerlist'] =  segmentmarkerlist
141    mesh_dict['triangleneighborlist'] = neighborlist
142    mesh_dict['qaz'] = 1 #debugging
143    ##print "r_test", r_test
144   
145    return mesh_dict
146
147def add_area_tag(regions):
148    """
149    So, what is the format?
150    A list with
151    [x,y,region_tag,area] OR [x,y,region_tag]
152    if it's [x,y,region_tag], add a 4th element, value of 0.0.
153    """
154    if isinstance(regions, ListType):
155        for i, region in enumerate(regions):
156            if len(region) == 3:
157                if isinstance(region, TupleType):
158                    #FIXME: How do you convert a tuple to a list?
159                    # I can do it a stupid way..
160                    tuple = region[:]
161                    regions[i] = []
162                    for j in tuple:
163                        regions[i].append(j)
164                    regions[i].append(0.0)
165                else:
166                    regions[i].append(0.0)
167                   
168                # let ensure numeric catch this..   
169                #len(region) <= 2:
170                #msg = 'ERROR: Inconsistent regions array.'
171                #raise msg
172            #elif
173    return regions
174
175if __name__ == "__main__":
176    pass 
Note: See TracBrowser for help on using the repository browser.