Skip to content

plotting extruded meshes - #5368

Open
danshapero wants to merge 12 commits into
releasefrom
drs/plotting
Open

plotting extruded meshes#5368
danshapero wants to merge 12 commits into
releasefrom
drs/plotting

Conversation

@danshapero

Copy link
Copy Markdown
Contributor

Fixes several plotting issues.

  • We can plot extruded meshes.
  • We can plot CircleManifoldMesh or an interval mesh lifted into 2D or 3D.
  • Remove some warnings related to missing subdomains.
  • Remove call to a deprecated matplotlib API for looking up color maps.
  • General cleanup of dead code.

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.

danshapero and others added 10 commits August 19, 2026 11:37
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`.
@danshapero danshapero changed the title Drs/plotting plotting extruded meshes Aug 20, 2026
@connorjward

Copy link
Copy Markdown
Contributor

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.

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 connorjward left a comment

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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.

Comment thread tests/firedrake/output/test_plotting.py Outdated


@pytest.mark.skipplot
def test_triplot_kwargs_not_mutated():

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Seems a bit of an odd thing to test

Comment thread firedrake/pyplot/mpl.py
Comment on lines +178 to +180
if mesh.extruded and mesh.variable_layers:
raise NotImplementedError("Visualizing variable layer extruded meshes "
"not implemented yet!")

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Suggested change
if mesh.extruded and mesh.variable_layers:
raise NotImplementedError("Visualizing variable layer extruded meshes "
"not implemented yet!")

variable layer extrusion is deprecated

Comment thread firedrake/pyplot/mpl.py Outdated
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>`

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Comment thread firedrake/pyplot/mpl.py
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

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

what about hexes or prisms?

Copy link
Copy Markdown
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

None yet

Projects

None yet

Development

Successfully merging this pull request may close these issues.

2 participants