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 allow_volume_steiner=True, quality_meshing=True,
139 generate_edges=None, generate_faces=False, min_angle=None):
140 """Triangulate the domain given in `mesh_info'."""
141 opts = "pzj"
142 if quality_meshing:
143 if min_angle is not None:
144 opts += "q%f" % min_angle
145 else:
146 opts += "q"
147
148 if verbose:
149 opts += "VV"
150 else:
151 opts += "Q"
152
153 if attributes:
154 opts += "A"
155
156 if volume_constraints:
157 opts += "a"
158 if max_volume:
159 raise ValueError, "cannot specify both volume_constraints and max_area"
160 elif max_volume:
161 opts += "a%s" % repr(max_volume)
162
163 if refinement_func is not None:
164 opts += "u"
165
166 if generate_edges is not None:
167 from warnings import warn
168 warn("generate_edges is deprecated--use generate_faces instead")
169 generate_faces = generate_edges
170
171 if generate_faces:
172 opts += "e"
173
174 if not allow_volume_steiner:
175 opts += "YY"
176 if allow_boundary_steiner:
177 raise ValueError("cannot allow boundary Steiner points when volume "
178 "Steiner points are forbidden")
179 else:
180 if not allow_boundary_steiner:
181 opts += "Y"
182
183
184 try:
185 import locale
186 except ImportErorr:
187 have_locale = False
188 else:
189 have_locale = True
190 prev_num_locale = locale.getlocale(locale.LC_NUMERIC)
191 locale.setlocale(locale.LC_NUMERIC, "C")
192
193 try:
194 mesh = MeshInfo()
195 internals.triangulate(opts, mesh_info, mesh, MeshInfo(), refinement_func)
196 finally:
197
198 if have_locale:
199 locale.setlocale(locale.LC_NUMERIC, prev_num_locale)
200
201 return mesh
202
203
204
205
206 -def refine(input_p, verbose=False, refinement_func=None):
207 opts = "razj"
208 if len(input_p.faces) != 0:
209 opts += "p"
210 if verbose:
211 opts += "VV"
212 else:
213 opts += "Q"
214 if refinement_func is not None:
215 opts += "u"
216 output_p = MeshInfo()
217 internals.triangulate(opts, input_p, output_p, MeshInfo(), refinement_func)
218 return output_p
219
220
221
222
224 gp_file = open(filename, "w")
225
226 if facets:
227 segments = out_p.facets
228 else:
229 segments = out_p.elements
230
231 for points in segments:
232 for pt in points:
233 gp_file.write("%f %f\n" % tuple(out_p.points[pt]))
234 gp_file.write("%f %f\n\n" % tuple(out_p.points[points[0]]))
235