Shortcuts on this page
r m x toggle line displays
j k next/prev highlighted chunk
0 (zero) top of page
1 (one) first highlighted chunk
1"""
2Mesh refinement for triangular grids.
3"""
5import numpy as np
7from matplotlib import _api
8from matplotlib.tri._triangulation import Triangulation
9import matplotlib.tri._triinterpolate
12class TriRefiner:
13 """
14 Abstract base class for classes implementing mesh refinement.
16 A TriRefiner encapsulates a Triangulation object and provides tools for
17 mesh refinement and interpolation.
19 Derived classes must implement:
21 - ``refine_triangulation(return_tri_index=False, **kwargs)`` , where
22 the optional keyword arguments *kwargs* are defined in each
23 TriRefiner concrete implementation, and which returns:
25 - a refined triangulation,
26 - optionally (depending on *return_tri_index*), for each
27 point of the refined triangulation: the index of
28 the initial triangulation triangle to which it belongs.
30 - ``refine_field(z, triinterpolator=None, **kwargs)``, where:
32 - *z* array of field values (to refine) defined at the base
33 triangulation nodes,
34 - *triinterpolator* is an optional `~matplotlib.tri.TriInterpolator`,
35 - the other optional keyword arguments *kwargs* are defined in
36 each TriRefiner concrete implementation;
38 and which returns (as a tuple) a refined triangular mesh and the
39 interpolated values of the field at the refined triangulation nodes.
40 """
42 def __init__(self, triangulation):
43 _api.check_isinstance(Triangulation, triangulation=triangulation)
44 self._triangulation = triangulation
47class UniformTriRefiner(TriRefiner):
48 """
49 Uniform mesh refinement by recursive subdivisions.
51 Parameters
52 ----------
53 triangulation : `~matplotlib.tri.Triangulation`
54 The encapsulated triangulation (to be refined)
55 """
56# See Also
57# --------
58# :class:`~matplotlib.tri.CubicTriInterpolator` and
59# :class:`~matplotlib.tri.TriAnalyzer`.
60# """
61 def __init__(self, triangulation):
62 super().__init__(triangulation)
64 def refine_triangulation(self, return_tri_index=False, subdiv=3):
65 """
66 Compute a uniformly refined triangulation *refi_triangulation* of
67 the encapsulated :attr:`triangulation`.
69 This function refines the encapsulated triangulation by splitting each
70 father triangle into 4 child sub-triangles built on the edges midside
71 nodes, recursing *subdiv* times. In the end, each triangle is hence
72 divided into ``4**subdiv`` child triangles.
74 Parameters
75 ----------
76 return_tri_index : bool, default: False
77 Whether an index table indicating the father triangle index of each
78 point is returned.
79 subdiv : int, default: 3
80 Recursion level for the subdivision.
81 Each triangle is divided into ``4**subdiv`` child triangles;
82 hence, the default results in 64 refined subtriangles for each
83 triangle of the initial triangulation.
85 Returns
86 -------
87 refi_triangulation : `~matplotlib.tri.Triangulation`
88 The refined triangulation.
89 found_index : int array
90 Index of the initial triangulation containing triangle, for each
91 point of *refi_triangulation*.
92 Returned only if *return_tri_index* is set to True.
93 """
94 refi_triangulation = self._triangulation
95 ntri = refi_triangulation.triangles.shape[0]
97 # Computes the triangulation ancestors numbers in the reference
98 # triangulation.
99 ancestors = np.arange(ntri, dtype=np.int32)
100 for _ in range(subdiv):
101 refi_triangulation, ancestors = self._refine_triangulation_once(
102 refi_triangulation, ancestors)
103 refi_npts = refi_triangulation.x.shape[0]
104 refi_triangles = refi_triangulation.triangles
106 # Now we compute found_index table if needed
107 if return_tri_index:
108 # We have to initialize found_index with -1 because some nodes
109 # may very well belong to no triangle at all, e.g., in case of
110 # Delaunay Triangulation with DuplicatePointWarning.
111 found_index = np.full(refi_npts, -1, dtype=np.int32)
112 tri_mask = self._triangulation.mask
113 if tri_mask is None:
114 found_index[refi_triangles] = np.repeat(ancestors,
115 3).reshape(-1, 3)
116 else:
117 # There is a subtlety here: we want to avoid whenever possible
118 # that refined points container is a masked triangle (which
119 # would result in artifacts in plots).
120 # So we impose the numbering from masked ancestors first,
121 # then overwrite it with unmasked ancestor numbers.
122 ancestor_mask = tri_mask[ancestors]
123 found_index[refi_triangles[ancestor_mask, :]
124 ] = np.repeat(ancestors[ancestor_mask],
125 3).reshape(-1, 3)
126 found_index[refi_triangles[~ancestor_mask, :]
127 ] = np.repeat(ancestors[~ancestor_mask],
128 3).reshape(-1, 3)
129 return refi_triangulation, found_index
130 else:
131 return refi_triangulation
133 def refine_field(self, z, triinterpolator=None, subdiv=3):
134 """
135 Refine a field defined on the encapsulated triangulation.
137 Parameters
138 ----------
139 z : (npoints,) array-like
140 Values of the field to refine, defined at the nodes of the
141 encapsulated triangulation. (``n_points`` is the number of points
142 in the initial triangulation)
143 triinterpolator : `~matplotlib.tri.TriInterpolator`, optional
144 Interpolator used for field interpolation. If not specified,
145 a `~matplotlib.tri.CubicTriInterpolator` will be used.
146 subdiv : int, default: 3
147 Recursion level for the subdivision.
148 Each triangle is divided into ``4**subdiv`` child triangles.
150 Returns
151 -------
152 refi_tri : `~matplotlib.tri.Triangulation`
153 The returned refined triangulation.
154 refi_z : 1D array of length: *refi_tri* node count.
155 The returned interpolated field (at *refi_tri* nodes).
156 """
157 if triinterpolator is None:
158 interp = matplotlib.tri.CubicTriInterpolator(
159 self._triangulation, z)
160 else:
161 _api.check_isinstance(matplotlib.tri.TriInterpolator,
162 triinterpolator=triinterpolator)
163 interp = triinterpolator
165 refi_tri, found_index = self.refine_triangulation(
166 subdiv=subdiv, return_tri_index=True)
167 refi_z = interp._interpolate_multikeys(
168 refi_tri.x, refi_tri.y, tri_index=found_index)[0]
169 return refi_tri, refi_z
171 @staticmethod
172 def _refine_triangulation_once(triangulation, ancestors=None):
173 """
174 Refine a `.Triangulation` by splitting each triangle into 4
175 child-masked_triangles built on the edges midside nodes.
177 Masked triangles, if present, are also split, but their children
178 returned masked.
180 If *ancestors* is not provided, returns only a new triangulation:
181 child_triangulation.
183 If the array-like key table *ancestor* is given, it shall be of shape
184 (ntri,) where ntri is the number of *triangulation* masked_triangles.
185 In this case, the function returns
186 (child_triangulation, child_ancestors)
187 child_ancestors is defined so that the 4 child masked_triangles share
188 the same index as their father: child_ancestors.shape = (4 * ntri,).
189 """
191 x = triangulation.x
192 y = triangulation.y
194 # According to tri.triangulation doc:
195 # neighbors[i, j] is the triangle that is the neighbor
196 # to the edge from point index masked_triangles[i, j] to point
197 # index masked_triangles[i, (j+1)%3].
198 neighbors = triangulation.neighbors
199 triangles = triangulation.triangles
200 npts = np.shape(x)[0]
201 ntri = np.shape(triangles)[0]
202 if ancestors is not None:
203 ancestors = np.asarray(ancestors)
204 if np.shape(ancestors) != (ntri,):
205 raise ValueError(
206 "Incompatible shapes provide for "
207 "triangulation.masked_triangles and ancestors: "
208 f"{np.shape(triangles)} and {np.shape(ancestors)}")
210 # Initiating tables refi_x and refi_y of the refined triangulation
211 # points
212 # hint: each apex is shared by 2 masked_triangles except the borders.
213 borders = np.sum(neighbors == -1)
214 added_pts = (3*ntri + borders) // 2
215 refi_npts = npts + added_pts
216 refi_x = np.zeros(refi_npts)
217 refi_y = np.zeros(refi_npts)
219 # First part of refi_x, refi_y is just the initial points
220 refi_x[:npts] = x
221 refi_y[:npts] = y
223 # Second part contains the edge midside nodes.
224 # Each edge belongs to 1 triangle (if border edge) or is shared by 2
225 # masked_triangles (interior edge).
226 # We first build 2 * ntri arrays of edge starting nodes (edge_elems,
227 # edge_apexes); we then extract only the masters to avoid overlaps.
228 # The so-called 'master' is the triangle with biggest index
229 # The 'slave' is the triangle with lower index
230 # (can be -1 if border edge)
231 # For slave and master we will identify the apex pointing to the edge
232 # start
233 edge_elems = np.tile(np.arange(ntri, dtype=np.int32), 3)
234 edge_apexes = np.repeat(np.arange(3, dtype=np.int32), ntri)
235 edge_neighbors = neighbors[edge_elems, edge_apexes]
236 mask_masters = (edge_elems > edge_neighbors)
238 # Identifying the "masters" and adding to refi_x, refi_y vec
239 masters = edge_elems[mask_masters]
240 apex_masters = edge_apexes[mask_masters]
241 x_add = (x[triangles[masters, apex_masters]] +
242 x[triangles[masters, (apex_masters+1) % 3]]) * 0.5
243 y_add = (y[triangles[masters, apex_masters]] +
244 y[triangles[masters, (apex_masters+1) % 3]]) * 0.5
245 refi_x[npts:] = x_add
246 refi_y[npts:] = y_add
248 # Building the new masked_triangles; each old masked_triangles hosts
249 # 4 new masked_triangles
250 # there are 6 pts to identify per 'old' triangle, 3 new_pt_corner and
251 # 3 new_pt_midside
252 new_pt_corner = triangles
254 # What is the index in refi_x, refi_y of point at middle of apex iapex
255 # of elem ielem ?
256 # If ielem is the apex master: simple count, given the way refi_x was
257 # built.
258 # If ielem is the apex slave: yet we do not know; but we will soon
259 # using the neighbors table.
260 new_pt_midside = np.empty([ntri, 3], dtype=np.int32)
261 cum_sum = npts
262 for imid in range(3):
263 mask_st_loc = (imid == apex_masters)
264 n_masters_loc = np.sum(mask_st_loc)
265 elem_masters_loc = masters[mask_st_loc]
266 new_pt_midside[:, imid][elem_masters_loc] = np.arange(
267 n_masters_loc, dtype=np.int32) + cum_sum
268 cum_sum += n_masters_loc
270 # Now dealing with slave elems.
271 # for each slave element we identify the master and then the inode
272 # once slave_masters is identified, slave_masters_apex is such that:
273 # neighbors[slaves_masters, slave_masters_apex] == slaves
274 mask_slaves = np.logical_not(mask_masters)
275 slaves = edge_elems[mask_slaves]
276 slaves_masters = edge_neighbors[mask_slaves]
277 diff_table = np.abs(neighbors[slaves_masters, :] -
278 np.outer(slaves, np.ones(3, dtype=np.int32)))
279 slave_masters_apex = np.argmin(diff_table, axis=1)
280 slaves_apex = edge_apexes[mask_slaves]
281 new_pt_midside[slaves, slaves_apex] = new_pt_midside[
282 slaves_masters, slave_masters_apex]
284 # Builds the 4 child masked_triangles
285 child_triangles = np.empty([ntri*4, 3], dtype=np.int32)
286 child_triangles[0::4, :] = np.vstack([
287 new_pt_corner[:, 0], new_pt_midside[:, 0],
288 new_pt_midside[:, 2]]).T
289 child_triangles[1::4, :] = np.vstack([
290 new_pt_corner[:, 1], new_pt_midside[:, 1],
291 new_pt_midside[:, 0]]).T
292 child_triangles[2::4, :] = np.vstack([
293 new_pt_corner[:, 2], new_pt_midside[:, 2],
294 new_pt_midside[:, 1]]).T
295 child_triangles[3::4, :] = np.vstack([
296 new_pt_midside[:, 0], new_pt_midside[:, 1],
297 new_pt_midside[:, 2]]).T
298 child_triangulation = Triangulation(refi_x, refi_y, child_triangles)
300 # Builds the child mask
301 if triangulation.mask is not None:
302 child_triangulation.set_mask(np.repeat(triangulation.mask, 4))
304 if ancestors is None:
305 return child_triangulation
306 else:
307 return child_triangulation, np.repeat(ancestors, 4)