source: inundation/pyvolution/pmesh2domain.py @ 1869

Last change on this file since 1869 was 1557, checked in by steve, 19 years ago
File size: 5.7 KB
Line 
1"""Class pmesh2domain - Converting .tsh files to doamains
2
3
4   Copyright 2004
5   Ole Nielsen, Stephen Roberts, Duncan Gray, Christopher Zoppou
6   Geoscience Australia
7"""
8
9
10def pmesh_to_domain_instance(fileName, DomainClass, setting_function = None):
11    """
12    """
13    import sys
14    from domain import Domain
15
16    vertex_coordinates, vertices, tag_dict, vertex_quantity_dict \
17                        ,tagged_elements_dict, geo_reference = \
18                pmesh_to_domain(fileName, setting_function=setting_function)
19
20
21    assert issubclass(DomainClass, Domain),"DomainClass is not a subclass of Domain."
22
23
24    domain = DomainClass(coordinates = vertex_coordinates,
25                         vertices = vertices,
26                         boundary = tag_dict,
27                         tagged_elements = tagged_elements_dict,
28                         geo_reference = geo_reference )
29
30
31
32    #FIXME (Ole): Is this really the right place to apply the a default
33    #value specific to the shallow water wave equation?
34    #The 'assert' above indicates that any subclass of Domain is acceptable.
35    #Suggestion - module shallow_water.py will eventually take care of this
36    #(when I get around to it) so it should be removed from here.
37
38    # This doesn't work on the domain instance.
39    # This is still needed so -ve elevations don't cuase 'lakes'
40    # The fixme we discussed was to only create a quantity when its values
41    #are set.
42    # I think that's the way to go still
43
44    # set the water stage to be the elevation
45    if vertex_quantity_dict.has_key('elevation') and not vertex_quantity_dict.has_key('stage'):
46        vertex_quantity_dict['stage'] = vertex_quantity_dict['elevation']
47
48    domain.set_quantity_vertices_dict(vertex_quantity_dict)
49    #print "vertex_quantity_dict",vertex_quantity_dict
50    return domain
51
52
53
54def pmesh_to_domain(fileName, setting_function = None):
55    """
56    convert a pmesh dictionary to a list of Volumes.
57    Also, return a list of triangles which have boundary tags
58    mesh_dict structure;
59    generated point list: [(x1,y1),(x2,y2),...] (Tuples of doubles)
60    generated point attribute list:[(P1att1,P1attt2, ...),(P2att1,P2attt2,...),...]
61    generated segment list: [(point1,point2),(p3,p4),...] (Tuples of integers)
62    generated segment tag list: [S1Tag, S2Tag, ...] (list of ints)
63    triangle list:  [(point1,point2, point3),(p5,p4, p1),...] (Tuples of integers)
64    triangle neighbor list: [(triangle1,triangle2, triangle3),(t5,t4, t1),...] (Tuples of integers) -1 means there's no triangle neighbor
65    triangle attribute list: [T1att, T2att, ...] (list of strings)
66    """
67
68    from Numeric import transpose
69    from load_mesh.loadASCII import import_mesh_file
70
71    mesh_dict = import_mesh_file(fileName)
72    #print "mesh_dict",mesh_dict
73    vertex_coordinates = mesh_dict['vertices']
74    volumes = mesh_dict['triangles']
75    vertex_quantity_dict = {}
76    point_atts = transpose(mesh_dict['vertex_attributes'])
77    point_titles  = mesh_dict['vertex_attribute_titles']
78    geo_reference  = mesh_dict['geo_reference']
79    for quantity, value_vector in map (None, point_titles, point_atts):
80        vertex_quantity_dict[quantity] = value_vector
81    tag_dict = pmesh_dict_to_tag_dict(mesh_dict)
82    tagged_elements_dict = build_tagged_elements_dictionary(mesh_dict)
83    return vertex_coordinates, volumes, tag_dict, vertex_quantity_dict, tagged_elements_dict, geo_reference
84
85
86
87def build_tagged_elements_dictionary(mesh_dict):
88    """Build the dictionary of element tags.
89    tagged_elements is a dictionary of element arrays,
90    keyed by tag:
91    { (tag): [e1, e2, e3..] }
92    """
93    tri_atts = mesh_dict['triangle_tags']
94    #print "tri_atts", tri_atts
95    tagged_elements = {}
96    for tri_att_index in range(len(tri_atts)):
97        tagged_elements.setdefault(tri_atts[tri_att_index],[]).append(tri_att_index)
98    #print "DSG pm2do tagged_elements", tagged_elements
99    return tagged_elements
100
101#FIXME: The issue is whether this format should be stored in the tsh file
102#instead of having to be created here?
103
104#This information is pyvolution focused, not mesh generation focused.
105#This is an appropriate place for the info to be created. -DSG
106
107#FIXME: Another issue is that the tsh file stores consecutive
108#indices explicitly. This is really redundant.
109#Suggest looking at obj and our own sww format and also consider
110#using netCDF.
111
112# It is redundant.  It was the format originally decided on.
113# I'm happy for it to change. -DSG
114
115def pmesh_dict_to_tag_dict(mesh_dict):
116    """ Convert the pmesh dictionary (mesh_dict) description of boundary tags
117    to a dictionary of tags, indexed with volume id and face number.
118    """
119    triangles = mesh_dict['triangles']
120    sides = calc_sides(triangles)
121    tag_dict = {}
122    for seg, tag in map(None,mesh_dict['segments'],
123                           mesh_dict['segment_tags']):
124        v1 = seg[0]
125        v2 = seg[1]
126        for key in [(v1,v2),(v2,v1)]:
127            if sides.has_key(key) and tag <> "":
128                #"" represents null.  Don't put these into the dictionary
129                #this creates a dict of lists of faces, indexed by tag
130                #tagged_edges.setdefault(tag,[]).append(sides[key])
131                tag_dict[sides[key]] = tag
132
133    return tag_dict
134
135
136def calc_sides(triangles):
137    #Build dictionary mapping from sides (2-tuple of points)
138    #to left hand side neighbouring triangle
139    sides = {}
140    for id, triangle in enumerate(triangles):
141        a = triangle[0]
142        b = triangle[1]
143        c = triangle[2]
144        sides[a,b] = (id, 2) #(id, face)
145        sides[b,c] = (id, 0) #(id, face)
146        sides[c,a] = (id, 1) #(id, face)
147    return sides
Note: See TracBrowser for help on using the repository browser.