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