plotting extruded meshes - #5368
Conversation
triplot rejects extruded meshes outright, silently mangles the quadrilateral facets of hexahedral meshes into triangles, and pops entries out of the caller's boundary_kw. FunctionPlotter assumes the vertical node offset of an extruded function space is 1, so it samples the wrong data for anything but CG1. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
triplot assumed that a facet has as many vertices as the mesh has topological dimensions, which is only true for simplices. The quadrilateral facets of a hexahedron were reshaped into triangles spanning several facets at once, so the boundary of a hexahedral mesh was drawn as a mess of overlapping slivers. Take the nodes of each facet from the topology of the reference cell instead, which gives the right answer for any cell. This also replaces the Python loop over facets with an indexing operation and hoists it out of the loop over markers, where it was being recomputed for every marker. While here, copy boundary_kw before popping the colours out of it, so that a dict of keyword arguments can be reused across calls. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
FunctionPlotter walked up a column of an extruded mesh by adding the layer number to the node numbers of the bottom cell. That is only right when consecutive layers are one node apart, which holds for CG1 but not for CG2 (offset 2), DQ1 (offset 4) or the coordinates of a periodic extrusion (offset 2). Every field that tricontourf, tricontour, tripcolor and trisurf drew on such a space was sampled from the wrong degrees of freedom. Use the offset that the cell node map already carries. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
An extruded mesh has two kinds of facet. The vertical facets are the extrusions of the facets of the base mesh and carry its markers, so they are drawn from the extruded facet node maps. The horizontal facets are not facets of the base mesh at all; they bound the bottom and top of each column of cells, and are drawn from the cell node map of the lowest and highest layer under the markers "bottom" and "top", matching the subdomain ids of the ds_b and ds_t measures. A periodic extrusion identifies the bottom with the top, so neither is a boundary and neither is drawn. The test that the coordinates are piecewise linear now compares the dimension of the coordinate space against the number of vertices of the cell. Comparing the degree against 1 fails on a tensor product element, whose degree is a tuple. Variable layer extrusion is still unsupported, but now says so instead of failing obscurely. It is deprecated for removal in 2026.10.0. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
triplot reached for matplotlib.cm.get_cmap first and only fell back to matplotlib.colormaps if it was missing, so every call raised a MatplotlibDeprecationWarning. The fallback is the newer of the two: matplotlib.colormaps has been available since 3.5, well before matplotlib.cm.get_cmap was deprecated in 3.7, so nothing needs the deprecated call and it can go. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
_entity_dimensions had one caller and _cycle had one caller, so neither earned a name of its own. Fold them into the call site and into _entity_node_list respectively, keeping what they documented. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
Repeating the first vertex of a cell was necessary when the interior was drawn with a LineCollection, which leaves its polylines open: without it every triangle came out missing an edge. The interior became a PolyCollection ten days later in 808ba85, and both it and Poly3DCollection close their polygons themselves, so the repeated vertex has been dead since. Dropping it leaves the 2D output pixel for pixel identical. In 3D it changes the depth sorting, which averages the z coordinates of a polygon's vertices: the repeated vertex weighted that average towards vertex 0, so the polygons are now sorted on their true centroid instead. The renderings are indistinguishable. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
The facets of a 1D mesh are points, which no line collection can show, so an interval embedded into 2D or 3D lost its two endpoints. Draw them with a PathCollection of markers instead, promoted to a Path3DCollection when the mesh lives in space. The markers are sized in points like those of a scatter plot, so only their offsets are in data coordinates and they need the axes to build. Choosing the collection from the geometric and topological dimensions cannot express that, since the types were looked up before the axes existed. Choose it from the shape of the vertex array instead: entities with one vertex are markers, those with two are segments, and those with more are polygons, with the trailing axis giving the dimension of the space they are drawn in. Two entities that agree on both are always drawn the same way, so the shape is all the dispatch needs, and it happens at the point of use where the axes are in hand. The rendering of every mesh that already plotted is unchanged, pixel for pixel. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
Every call to triplot was producing warnings about empty subsets because of how the querying for marked facets worked. This patch changes triplot to only color marked exterior facets, which gets rid of the warning and simplifies the code. But it won't plot any marked interior facets. Fixing that requires some changes to `mesh.py`.
Yes sorry about those warnings. I have looked into them getting raised in this file before and couldn't immediately figure out a nice fix. It also intersects with pyop3 changes and so I didn't want to fix something that will immediately need changing again. |
connorjward
left a comment
There was a problem hiding this comment.
The broken links are not your concern, but the docs are failing for legitimate reasons.
Also should we amend the plotting docs anywhere? Possibly not.
|
|
||
|
|
||
| @pytest.mark.skipplot | ||
| def test_triplot_kwargs_not_mutated(): |
There was a problem hiding this comment.
Seems a bit of an odd thing to test
| if mesh.extruded and mesh.variable_layers: | ||
| raise NotImplementedError("Visualizing variable layer extruded meshes " | ||
| "not implemented yet!") |
There was a problem hiding this comment.
| if mesh.extruded and mesh.variable_layers: | |
| raise NotImplementedError("Visualizing variable layer extruded meshes " | |
| "not implemented yet!") |
variable layer extrusion is deprecated
| that a mesh of any dimension is drawn from. | ||
|
|
||
| :arg vertices: array of shape ``(num_entities, num_vertices, gdim)`` | ||
| :return: the matplotlib :class:`Collection <matplotlib.collections.Collection>` |
There was a problem hiding this comment.
Our style guide now wants numpydoc format: https://github.com/firedrakeproject/firedrake/wiki/Firedrake-Coding-Guide
| The vertices of each entity come out in a cycle around its perimeter. FInAT | ||
| numbers the vertices of a quadrilateral lexicographically, putting 1 and 2 | ||
| diagonally opposite each other; exchanging the last two makes the polygon | ||
| simple. Entities with any other number of vertices are simplices, which are |
There was a problem hiding this comment.
what about hexes or prisms?
There was a problem hiding this comment.
We only plot 3D meshes by drawing their faces, not their volumes, so this code is only called for triangles and quads. If we did some kind of fancy volume rendering that we'd have to worry about hexes or prisms.
Fixes several plotting issues.
I think all of this is straightforward except for how to remove the warnings related to missing subdomains is in the last commit. Before, we were getting all unique markers found throughout the whole mesh and querying which markers were present in either the interior or exterior facets. This was there in order to draw the boundary facets with colors. But if you ask "give me all the interior facets have marker X" when there are no interior facets that have marker X, or any marker at all, you get this warning about an empty subdomain. The most common use case is that there are no marked interior facets at all, so I opted to silence the warning in that case. Fixing this in general is going to require some changes in mesh.py which should go in a separate PR.
The code is all vibed but I checked the outputs myself.