1 from __future__ import division
2 from meshpy.common import MeshInfoBase, dump_array
3 import meshpy._triangle as internals
4
5
6
7
8 -class MeshInfo(internals.MeshInfo, MeshInfoBase):
9 _constituents = [
10 "points", "point_attributes", "point_markers",
11 "elements", "element_attributes", "element_volumes",
12 "neighbors",
13 "facets", "facet_markers",
14 "holes",
15 "regions",
16 "faces", "face_markers",
17 "normals",
18 ]
19
21 def dump_array(array):
22 try:
23 return [
24 [array[i,j] for j in range(array.unit)]
25 for i in range(len(array))]
26 except RuntimeError:
27
28 return None
29
30 return self.number_of_point_attributes, \
31 self.number_of_element_attributes, \
32 [(name, dump_array(getattr(self, name))) for name in self._constituents]
33
34 - def __setstate__(self, (p_attr_count, e_attr_count, state)):
35 self.number_of_point_attributes = p_attr_count
36 self.number_of_element_attributes = e_attr_count
37 for name, array in state:
38 if name not in self._constituents:
39 raise RuntimeError, "Unknown constituent during unpickling"
40
41 dest_array = getattr(self, name)
42
43 if array is None:
44 dest_array.deallocate()
45 else:
46 if len(dest_array) != len(array):
47 dest_array.resize(len(array))
48 if not dest_array.allocated and len(array)>0:
49 dest_array.setup()
50
51 for i,tup in enumerate(array):
52 for j,v in enumerate(tup):
53 dest_array[i, j] = v
54
56 self.facets.resize(len(facets))
57
58 for i, facet in enumerate(facets):
59 self.facets[i] = facet
60
61 if facet_markers is not None:
62 self.facet_markers.setup()
63 for i, mark in enumerate(facet_markers):
64 self.facet_markers[i] = mark
65
69
70
71
72
74 """Return a new facets array in which the original facets are
75 each subdivided into C{subdivisions} subfacets.
76
77 This routine is useful if you have to prohibit the insertion of Steiner
78 points on the boundary of your triangulation to allow the mesh to conform
79 either to itself periodically or another given mesh. In this case, you may
80 use this routine to create the necessary resolution along the boundary
81 in a predefined way.
82
83 @arg subdivisions: Either an C{int}, indicating a uniform number of subdivisions
84 throughout, or a list of the same length as C{facets}, specifying a subdivision
85 count for each individual facet.
86 @arg points: A list of points referred to from the facets list.
87 @arg facets: The list of old facets, in the form C{[(p1, p2), (p3,p4), ...]}.
88 @arg facet_markers: Either C{None} or a list of facet markers of the same length
89 as C{facets}.
90 @return: The new tuple C{(new_points, new_facets)}.
91 (Or C{(new_points, new_facets, new_facet_markers)} if C{facet_markers} is not
92 C{None}.)
93 """
94
95 def intermediate_points(pa, pb, n):
96 for i in range(1, n):
97 tau = i/n
98 yield [pai*(1-tau) + tau*pbi for pai, pbi in zip(pa, pb)]
99
100 if isinstance(subdivisions, int):
101 from itertools import repeat
102 subdiv_it = repeat(subdivisions, len(facets))
103 else:
104 assert len(facets) == len(subdivisions)
105 subdiv_it = subdivisions.__iter__()
106
107 new_points = points[:]
108 new_facets = []
109
110 if facet_markers is not None:
111 assert len(facets) == len(facet_markers)
112 new_facet_markers = []
113
114 for facet_idx, ((pidx_a, pidx_b), subdiv) in enumerate(zip(facets, subdiv_it)):
115 facet_points = [pidx_a]
116 for p in intermediate_points(points[pidx_a], points[pidx_b], subdiv):
117 facet_points.append(len(new_points))
118 new_points.append(p)
119 facet_points.append(pidx_b)
120
121 for i, p1 in enumerate(facet_points[:-1]):
122 p2 = facet_points[i+1]
123 new_facets.append((p1, p2))
124
125 if facet_markers is not None:
126 new_facet_markers.append(facet_markers[facet_idx])
127
128 if facet_markers is not None:
129 return new_points, new_facets, new_facet_markers
130 else:
131 return new_points, new_facets
132
133
134
135
136 -def build(mesh_info, verbose=False, refinement_func=None, attributes=False,
137 volume_constraints=False, max_volume=None, allow_boundary_steiner=True,
138 generate_edges=None, generate_faces=False):
139 """Triangulate the domain given in `mesh_info'."""
140 opts = "pzqj"
141 if verbose:
142 opts += "VV"
143 else:
144 opts += "Q"
145
146 if attributes:
147 opts += "A"
148
149 if volume_constraints:
150 opts += "a"
151 if max_volume:
152 raise ValueError, "cannot specify both volume_constraints and max_area"
153 elif max_volume:
154 opts += "a%s" % repr(max_volume)
155
156 if refinement_func is not None:
157 opts += "u"
158
159 if generate_edges is not None:
160 from warnings import warn
161 warn("generate_edges is deprecated--use generate_faces instead")
162 generate_faces = generate_edges
163
164 if generate_faces:
165 opts += "e"
166
167 if not allow_boundary_steiner:
168 opts += "Y"
169
170
171 try:
172 import locale
173 except ImportErorr:
174 have_locale = False
175 else:
176 have_locale = True
177 prev_num_locale = locale.getlocale(locale.LC_NUMERIC)
178 locale.setlocale(locale.LC_NUMERIC, "C")
179
180 try:
181 mesh = MeshInfo()
182 internals.triangulate(opts, mesh_info, mesh, MeshInfo(), refinement_func)
183 finally:
184
185 if have_locale:
186 locale.setlocale(locale.LC_NUMERIC, prev_num_locale)
187
188 return mesh
189
190
191
192
193 -def refine(input_p, verbose=False, refinement_func=None):
194 opts = "razj"
195 if len(input_p.faces) != 0:
196 opts += "p"
197 if verbose:
198 opts += "VV"
199 else:
200 opts += "Q"
201 if refinement_func is not None:
202 opts += "u"
203 output_p = MeshInfo()
204 internals.triangulate(opts, input_p, output_p, MeshInfo(), refinement_func)
205 return output_p
206
207
208
209
211 gp_file = open(filename, "w")
212
213 if facets:
214 segments = out_p.facets
215 else:
216 segments = out_p.elements
217
218 for points in segments:
219 for pt in points:
220 gp_file.write("%f %f\n" % tuple(out_p.points[pt]))
221 gp_file.write("%f %f\n\n" % tuple(out_p.points[points[0]]))
222