From d19d6c8bb9a08b41b6c38757d3c95ab03d2d7dee Mon Sep 17 00:00:00 2001 From: jorgensd Date: Wed, 16 Sep 2026 22:49:39 +0200 Subject: [PATCH 1/5] Add parallel plotting in single window --- chapter1/complex_mode.py | 51 +++++++++++----- chapter1/fundamentals_code.py | 84 +++++++++++++++++++------ chapter1/membrane_code.py | 69 +++++++++++++++------ chapter1/nitsche.py | 30 ++++++--- chapter2/amr.py | 84 ++++++++++++++++--------- chapter2/diffusion_code.py | 63 ++++++++++++------- chapter2/hyperelasticity.py | 61 ++++++++++++------ chapter2/linearelasticity_code.py | 64 +++++++++++++------ chapter2/ns_code1.py | 34 ++++++++--- chapter3/component_bc.py | 31 +++++++--- chapter3/em.py | 95 +++++++++++++++++++---------- chapter3/multiple_dirichlet.py | 31 +++++++--- chapter3/neumann_dirichlet_code.py | 33 +++++++--- chapter3/robin_neumann_dirichlet.py | 31 +++++++--- chapter3/subdomains.py | 91 +++++++++++++++++++-------- chapter4/newton-solver.py | 26 ++++++-- 16 files changed, 627 insertions(+), 251 deletions(-) diff --git a/chapter1/complex_mode.py b/chapter1/complex_mode.py index 83e8a2a8..e25c5323 100644 --- a/chapter1/complex_mode.py +++ b/chapter1/complex_mode.py @@ -180,6 +180,10 @@ # Finally, we plot the real and imaginary solutions. # +# Each process builds the grid of the cells it owns, the pieces are gathered on one process and +# merged into a single grid, which is then drawn, see +# [Plotting in parallel](./fundamentals_code). + # + import pyvista @@ -188,19 +192,34 @@ grid = pyvista.UnstructuredGrid(pyvista_cells, cell_types, geometry) grid.point_data["u_real"] = uh.x.array.real grid.point_data["u_imag"] = uh.x.array.imag -_ = grid.set_active_scalars("u_real") - -p_real = pyvista.Plotter() -p_real.add_text("uh real", position="upper_edge", font_size=14, color="black") -p_real.add_mesh(grid, show_edges=True) -p_real.view_xy() -if not pyvista.OFF_SCREEN: - p_real.show() - -grid.set_active_scalars("u_imag") -p_imag = pyvista.Plotter() -p_imag.add_text("uh imag", position="upper_edge", font_size=14, color="black") -p_imag.add_mesh(grid, show_edges=True) -p_imag.view_xy() -if not pyvista.OFF_SCREEN: - p_imag.show() + +root = 0 +assert root < mesh.comm.size, f"Cannot gather on process {root} of {mesh.comm.size}" + +clim_real = [ + mesh.comm.reduce(uh.x.array.real.min(), op=MPI.MIN, root=root), + mesh.comm.reduce(uh.x.array.real.max(), op=MPI.MAX, root=root), +] +clim_imag = [ + mesh.comm.reduce(uh.x.array.imag.min(), op=MPI.MIN, root=root), + mesh.comm.reduce(uh.x.array.imag.max(), op=MPI.MAX, root=root), +] +pieces = mesh.comm.gather(grid, root=root) +# - + +if pieces is not None: + merged_grid = pyvista.merge(pieces) + + p_real = pyvista.Plotter() + p_real.add_text("uh real", position="upper_edge", font_size=14, color="black") + p_real.add_mesh(merged_grid, scalars="u_real", show_edges=True, clim=clim_real) + p_real.view_xy() + if not pyvista.OFF_SCREEN: + p_real.show() + + p_imag = pyvista.Plotter() + p_imag.add_text("uh imag", position="upper_edge", font_size=14, color="black") + p_imag.add_mesh(merged_grid, scalars="u_imag", show_edges=True, clim=clim_imag) + p_imag.view_xy() + if not pyvista.OFF_SCREEN: + p_imag.show() diff --git a/chapter1/fundamentals_code.py b/chapter1/fundamentals_code.py index 4707befa..018b26f5 100644 --- a/chapter1/fundamentals_code.py +++ b/chapter1/fundamentals_code.py @@ -317,6 +317,26 @@ import pyvista print(pyvista.global_theme.jupyter_backend) +# - + +# ### Plotting in parallel +# When the program runs on several processes, each process holds one piece of the mesh, and knows +# nothing about the rest of it, so no process can draw the whole figure on its own. +# We therefore let each process build the grid of its own piece, send the pieces to a single +# process with {py:meth}`gather`, and let that process join them into a +# single grid with {py:func}`pyvista.merge`. +# {py:func}`dolfinx.plot.vtk_mesh` uses only the cells that the process *owns*, so a cell that is +# shared between two processes is not drawn twice. +# `merge` welds the points that the pieces have in common back together, so the merged grid behaves +# like one built on a single process: filters that follow the field from one cell into the next, +# such as {py:meth}`streamlines`, do not stop at the partition +# boundaries. +# +# The process that collects the pieces is chosen with `root`, which has to be one of the processes +# we run on. +# `gather` returns the list of pieces on that process, and `None` on every other process, so the +# returned value tells each process whether it is the one that draws. +# The same code runs on a single process, where the list holds a single piece. # + from dolfinx import plot @@ -324,6 +344,11 @@ domain.topology.create_connectivity(tdim, tdim) topology, cell_types, geometry = plot.vtk_mesh(domain, tdim) grid = pyvista.UnstructuredGrid(topology, cell_types, geometry) + +root = 0 +assert root < domain.comm.size, f"Cannot gather on process {root} of {domain.comm.size}" + +pieces = domain.comm.gather(grid, root=root) # - # There are several backends that can be used with pyvista, and they have different benefits and drawbacks. @@ -334,13 +359,15 @@ # In the jupyter notebook environment, we use the default setting of `pyvista.OFF_SCREEN=False`, # which will render plots directly in the notebook. -plotter = pyvista.Plotter() -plotter.add_mesh(grid, show_edges=True) -plotter.view_xy() -if not pyvista.OFF_SCREEN: - plotter.show() -else: - figure = plotter.screenshot("fundamentals_mesh.png") +if pieces is not None: + merged_grid = pyvista.merge(pieces) + plotter = pyvista.Plotter() + plotter.add_mesh(merged_grid, show_edges=True) + plotter.view_xy() + if not pyvista.OFF_SCREEN: + plotter.show() + else: + figure = plotter.screenshot("fundamentals_mesh.png") # ## Plotting a function using pyvista # We want to plot the solution `uh`. @@ -352,23 +379,44 @@ u_topology, u_cell_types, u_geometry = plot.vtk_mesh(V) # Next, we create the {py:class}`pyvista.UnstructuredGrid` and add the dof-values to the mesh. +# The values are attached before the grids are gathered, so that each piece carries its own data +# into the merged grid. +# We also compute the range of `uh` over all processes with +# {py:meth}`reduce`, and pass it to the plotter as `clim`. +# Only the root process draws, so the result is only needed there, and `reduce` leaves it as `None` +# on the other processes. +# + u_grid = pyvista.UnstructuredGrid(u_topology, u_cell_types, u_geometry) u_grid.point_data["u"] = uh.x.array.real u_grid.set_active_scalars("u") -u_plotter = pyvista.Plotter() -u_plotter.add_mesh(u_grid, show_edges=True) -u_plotter.view_xy() -if not pyvista.OFF_SCREEN: - u_plotter.show() -# We can also warp the mesh by scalar to make use of the 3D plotting. +u_clim = [ + domain.comm.reduce(uh.x.array.real.min(), op=MPI.MIN, root=root), + domain.comm.reduce(uh.x.array.real.max(), op=MPI.MAX, root=root), +] +u_pieces = domain.comm.gather(u_grid, root=root) +# - -warped = u_grid.warp_by_scalar() -plotter2 = pyvista.Plotter() -plotter2.add_mesh(warped, show_edges=True, show_scalar_bar=True) -if not pyvista.OFF_SCREEN: - plotter2.show() +if u_pieces is not None: + merged_u_grid = pyvista.merge(u_pieces) + u_plotter = pyvista.Plotter() + u_plotter.add_mesh(merged_u_grid, show_edges=True, clim=u_clim) + u_plotter.view_xy() + if not pyvista.OFF_SCREEN: + u_plotter.show() + +# We can also warp the mesh by scalar to make use of the 3D plotting. +# As the grid has already been merged, filters such as +# {py:meth}`warp_by_scalar` are applied to the whole mesh. + +if u_pieces is not None: + plotter2 = pyvista.Plotter() + plotter2.add_mesh( + merged_u_grid.warp_by_scalar(), show_edges=True, show_scalar_bar=True, clim=u_clim + ) + if not pyvista.OFF_SCREEN: + plotter2.show() # ## External post-processing # For post-processing outside the python code, it is suggested to save the solution to file using either diff --git a/chapter1/membrane_code.py b/chapter1/membrane_code.py index 0d2584c5..0612710e 100644 --- a/chapter1/membrane_code.py +++ b/chapter1/membrane_code.py @@ -180,7 +180,10 @@ def on_boundary(x): from dolfinx.plot import vtk_mesh import pyvista -# Extract topology from mesh and create {py:class}`pyvista.UnstructuredGrid` +# Extract topology from mesh and create {py:class}`pyvista.UnstructuredGrid`. +# Each process builds the grid of the cells it owns, the pieces are gathered on one process and +# merged into a single grid, which is then drawn, see +# [Plotting in parallel](./fundamentals_code). topology, cell_types, x = vtk_mesh(V) grid = pyvista.UnstructuredGrid(topology, cell_types, x) @@ -189,29 +192,59 @@ def on_boundary(x): # + grid.point_data["u"] = uh.x.array -warped = grid.warp_by_scalar("u", factor=25) - -plotter = pyvista.Plotter() -plotter.add_mesh(warped, show_edges=True, show_scalar_bar=True, scalars="u") -if not pyvista.OFF_SCREEN: - plotter.show() -else: - plotter.screenshot("deflection.png") + +root = 0 +assert root < domain.comm.size, f"Cannot gather on process {root} of {domain.comm.size}" + +u_clim = [ + domain.comm.reduce(uh.x.array.real.min(), op=MPI.MIN, root=root), + domain.comm.reduce(uh.x.array.real.max(), op=MPI.MAX, root=root), +] +u_pieces = domain.comm.gather(grid, root=root) + +if u_pieces is not None: + merged_grid = pyvista.merge(u_pieces) + plotter = pyvista.Plotter() + plotter.add_mesh( + merged_grid.warp_by_scalar("u", factor=25), + show_edges=True, + show_scalar_bar=True, + scalars="u", + clim=u_clim, + ) + if not pyvista.OFF_SCREEN: + plotter.show() + else: + plotter.screenshot("deflection.png") # - # We next plot the load on the domain -load_plotter = pyvista.Plotter() +# + p_grid = pyvista.UnstructuredGrid(*vtk_mesh(Q)) p_grid.point_data["p"] = pressure.x.array.real -warped_p = p_grid.warp_by_scalar("p", factor=0.5) -warped_p.set_active_scalars("p") -load_plotter.add_mesh(warped_p, show_scalar_bar=True) -load_plotter.view_xy() -if not pyvista.OFF_SCREEN: - load_plotter.show() -else: - load_plotter.screenshot("load.png") + +p_clim = [ + domain.comm.reduce(pressure.x.array.real.min(), op=MPI.MIN, root=root), + domain.comm.reduce(pressure.x.array.real.max(), op=MPI.MAX, root=root), +] +p_pieces = domain.comm.gather(p_grid, root=root) + +if p_pieces is not None: + merged_p_grid = pyvista.merge(p_pieces) + load_plotter = pyvista.Plotter() + load_plotter.add_mesh( + merged_p_grid.warp_by_scalar("p", factor=0.5), + show_scalar_bar=True, + scalars="p", + clim=p_clim, + ) + load_plotter.view_xy() + if not pyvista.OFF_SCREEN: + load_plotter.show() + else: + load_plotter.screenshot("load.png") +# - # ## Making curve plots throughout the domain # Another way to compare the deflection and the load is to make a plot along the line $x=0$. diff --git a/chapter1/nitsche.py b/chapter1/nitsche.py index 7b3b69c6..cd507ca3 100644 --- a/chapter1/nitsche.py +++ b/chapter1/nitsche.py @@ -127,19 +127,35 @@ # we no longer fullfill the equation to machine precision at the mesh vertices. # We also plot the solution using {py:mod}`pyvista` +# Each process builds the grid of the cells it owns, the pieces are gathered on one process and +# merged into a single grid, which is then drawn, see +# [Plotting in parallel](./fundamentals_code). + # + import pyvista grid = pyvista.UnstructuredGrid(*plot.vtk_mesh(V)) grid.point_data["u"] = uh.x.array.real grid.set_active_scalars("u") -plotter = pyvista.Plotter() -plotter.add_mesh(grid, show_edges=True, show_scalar_bar=True) -plotter.view_xy() -if not pyvista.OFF_SCREEN: - plotter.show() -else: - figure = plotter.screenshot("nitsche.png") + +root = 0 +assert root < domain.comm.size, f"Cannot gather on process {root} of {domain.comm.size}" + +clim = [ + domain.comm.reduce(uh.x.array.real.min(), op=MPI.MIN, root=root), + domain.comm.reduce(uh.x.array.real.max(), op=MPI.MAX, root=root), +] +pieces = domain.comm.gather(grid, root=root) + +if pieces is not None: + grid = pyvista.merge(pieces) + plotter = pyvista.Plotter() + plotter.add_mesh(grid, show_edges=True, show_scalar_bar=True, clim=clim) + plotter.view_xy() + if not pyvista.OFF_SCREEN: + plotter.show() + else: + figure = plotter.screenshot("nitsche.png") # - # ```{bibliography} diff --git a/chapter2/amr.py b/chapter2/amr.py index 3e38edef..0ad3c892 100644 --- a/chapter2/amr.py +++ b/chapter2/amr.py @@ -71,17 +71,27 @@ # We use pyvista to visualize the mesh. +# Each process builds the grid of the cells it owns, the pieces are gathered on one process and +# merged into a single grid, which is then drawn, see +# [Plotting in parallel](../chapter1/fundamentals_code). + # + tags=["hide-input"] grid = pyvista.UnstructuredGrid(*dolfinx.plot.vtk_mesh(mesh)) -grid.cell_data["ct"] = ct.values - -plotter = pyvista.Plotter() -plotter.add_mesh( - grid, show_edges=True, scalars="ct", cmap="blues", show_scalar_bar=False -) -plotter.view_xy() -if not pyvista.OFF_SCREEN: - plotter.show() +grid.cell_data["ct"] = ct.values[: grid.n_cells] + +root = 0 +assert root < mesh.comm.size, f"Cannot gather on process {root} of {mesh.comm.size}" + +pieces = mesh.comm.gather(grid, root=root) +if pieces is not None: + merged_grid = pyvista.merge(pieces) + plotter = pyvista.Plotter() + plotter.add_mesh( + merged_grid, show_edges=True, scalars="ct", cmap="blues", show_scalar_bar=False + ) + plotter.view_xy() + if not pyvista.OFF_SCREEN: + plotter.show() # - # We have read in any cell and facet markers that have been defined in the NetGen model, @@ -96,15 +106,23 @@ # + tags=["hide-input"] curved_grid = pyvista.UnstructuredGrid(*dolfinx.plot.vtk_mesh(curved_mesh)) -curved_grid.cell_data["ct"] = ct.values -plotter = pyvista.Plotter() -plotter.add_mesh( - curved_grid, show_edges=False, scalars="ct", cmap="blues", show_scalar_bar=False -) -plotter.add_mesh(grid, style="wireframe", color="black") -plotter.view_xy() -if not pyvista.OFF_SCREEN: - plotter.show() +curved_grid.cell_data["ct"] = ct.values[: curved_grid.n_cells] + +curved_pieces = curved_mesh.comm.gather(curved_grid, root=root) +if curved_pieces is not None and pieces is not None: + merged_curved_grid = pyvista.merge(curved_pieces) + plotter = pyvista.Plotter() + plotter.add_mesh( + merged_curved_grid, + show_edges=False, + scalars="ct", + cmap="blues", + show_scalar_bar=False, + ) + plotter.add_mesh(merged_grid, style="wireframe", color="black") + plotter.view_xy() + if not pyvista.OFF_SCREEN: + plotter.show() # - # ## Solving the eigenvalue problem @@ -238,15 +256,19 @@ def mark_cells(uh_r: dolfinx.fem.Function, lam: float): # We will track the progress of the adaptive mesh refinement as a GIF. -plotter = pyvista.Plotter() -plotter.open_gif("amr.gif", fps=1) +plotter = None +if mesh.comm.rank == 0: + plotter = pyvista.Plotter() + plotter.open_gif("amr.gif", fps=1) # We make a convenience function to attach the relevant data to the plotter at a given # refinement step. +# The mesh is refined between the frames, so unlike the other animations in this tutorial, the +# pieces are gathered anew at each frame. # + tags=["hide-input"] -def write_frame(plotter: pyvista.Plotter, uh_r: dolfinx.fem.Function): +def write_frame(plotter: pyvista.Plotter | None, uh_r: dolfinx.fem.Function): # Scale uh_r to be consistent between refinement steps, as it can be multiplied by -1 uh_r_min = curved_mesh.comm.allreduce(uh_r.x.array.min(), op=MPI.MIN) uh_r_max = curved_mesh.comm.allreduce(uh_r.x.array.max(), op=MPI.MAX) @@ -261,16 +283,21 @@ def write_frame(plotter: pyvista.Plotter, uh_r: dolfinx.fem.Function): curved_grid = pyvista.UnstructuredGrid(*dolfinx.plot.vtk_mesh(uh_r.function_space)) curved_grid.point_data["u"] = uh_r.x.array curved_grid = curved_grid.tessellate() - curved_actor = plotter.add_mesh( - curved_grid, + + # The gathers are collective, so every process takes part in them, while only the process + # that received the pieces draws the frame + pieces = mesh.comm.gather(grid, root=root) + curved_pieces = mesh.comm.gather(curved_grid, root=root) + if plotter is None or pieces is None or curved_pieces is None: + return + plotter.clear() + plotter.add_mesh( + pyvista.merge(curved_pieces), show_edges=False, ) - - actor = plotter.add_mesh(grid, style="wireframe", color="black") + plotter.add_mesh(pyvista.merge(pieces), style="wireframe", color="black") plotter.view_xy() plotter.write_frame() - plotter.remove_actor(actor) - plotter.remove_actor(curved_actor) # - @@ -304,7 +331,8 @@ def write_frame(plotter: pyvista.Plotter, uh_r: dolfinx.fem.Function): if relative_error < termination_criteria: PETSc.Sys.Print(f"Converged in {i + 1} iterations.") break -plotter.close() +if plotter is not None: + plotter.close() # - # gif diff --git a/chapter2/diffusion_code.py b/chapter2/diffusion_code.py index 1593521a..183a9bb3 100644 --- a/chapter2/diffusion_code.py +++ b/chapter2/diffusion_code.py @@ -168,14 +168,18 @@ def initial_condition(x, a=5): # We would also like to visualize a colorbar reflecting the minimal and maximum # value of $u$ at each time step. -# + -grid = pyvista.UnstructuredGrid(*plot.vtk_mesh(V)) +# The mesh does not change during the simulation, so the pieces are gathered once, and only the +# values are sent to the process that writes the GIF at each frame. +# That process attaches the values to the pieces and merges them into a single grid before drawing, +# see [Plotting in parallel](../chapter1/fundamentals_code). +# We collect the pieces on the process `root`, which has to be one of the processes we run on. -plotter = pyvista.Plotter() -plotter.open_gif("u_time.gif", fps=10) +# + +root = 0 +assert root < domain.comm.size, f"Cannot gather on process {root} of {domain.comm.size}" -grid.point_data["uh"] = uh.x.array -warped = grid.warp_by_scalar("uh", factor=1) +grid = pyvista.UnstructuredGrid(*plot.vtk_mesh(V)) +pieces = domain.comm.gather(grid, root=root) viridis = mpl.colormaps.get_cmap("viridis").resampled(25) sargs = dict( @@ -188,15 +192,30 @@ def initial_condition(x, a=5): width=0.8, height=0.1, ) - -renderer = plotter.add_mesh( - warped, - show_edges=True, - lighting=False, - cmap=viridis, - scalar_bar_args=sargs, - clim=[0, max(uh.x.array)], -) +clim = [0, domain.comm.reduce(uh.x.array.max(), op=MPI.MAX, root=root)] + +plotter = None +if pieces is not None: + plotter = pyvista.Plotter() + plotter.open_gif("u_time.gif", fps=10) + + +def write_frame(values: list[np.ndarray]): + """Draw the merged solution into a single frame of the GIF.""" + assert plotter is not None and pieces is not None + for piece, piece_values in zip(pieces, values): + piece.point_data["uh"] = piece_values + merged_grid = pyvista.merge(pieces) + plotter.clear() + plotter.add_mesh( + merged_grid.warp_by_scalar("uh", factor=1), + show_edges=True, + lighting=False, + cmap=viridis, + scalar_bar_args=sargs, + clim=clim, + ) + plotter.write_frame() # - # (time-dep-assembly)= @@ -242,12 +261,14 @@ def initial_condition(x, a=5): # Write solution to file xdmf.write_function(uh, t) - # Update plot - new_warped = grid.warp_by_scalar("uh", factor=1) - warped.points[:, :] = new_warped.points - warped.point_data["uh"][:] = uh.x.array - plotter.write_frame() -plotter.close() + + # Update plot. The gather is collective, so every process takes part in it, while only the + # process that received the pieces draws the frame + values = domain.comm.gather(uh.x.array.copy(), root=root) + if values is not None: + write_frame(values) +if plotter is not None: + plotter.close() xdmf.close() # We {py:meth}`destroy` the PETSc objects to avoid memory leaks. diff --git a/chapter2/hyperelasticity.py b/chapter2/hyperelasticity.py index 17572923..2559a438 100644 --- a/chapter2/hyperelasticity.py +++ b/chapter2/hyperelasticity.py @@ -173,10 +173,14 @@ def right(x): # We create a function to plot the solution at each time step. -# + -plotter = pyvista.Plotter() -plotter.open_gif("deformation.gif", fps=3) +# Each process builds the grid of the cells it owns, the pieces are gathered on one process and +# merged into a single grid, which is then drawn, see +# [Plotting in parallel](../chapter1/fundamentals_code). +# +# The mesh does not change during the simulation, so the pieces are gathered once, and only the +# values are sent to the process that writes the GIF at each frame. +# + topology, cells, geometry = plot.vtk_mesh(u.function_space) function_grid = pyvista.UnstructuredGrid(topology, cells, geometry) @@ -185,13 +189,6 @@ def right(x): function_grid["u"] = values function_grid.set_active_vectors("u") -# Warp mesh by deformation -warped = function_grid.warp_by_vector("u", factor=1) -warped.set_active_vectors("u") - -# Add mesh to plotter and visualize -actor = plotter.add_mesh(warped, show_edges=True, lighting=False, clim=[0, 10]) - # Compute magnitude of displacement to visualize in GIF Vs = fem.functionspace(domain, ("Lagrange", 2)) magnitude = fem.Function(Vs) @@ -199,7 +196,31 @@ def right(x): ufl.sqrt(sum([u[i] ** 2 for i in range(len(u))])), Vs.element.interpolation_points ) magnitude.interpolate(us) -warped["mag"] = magnitude.x.array +function_grid["mag"] = magnitude.x.array + +root = 0 +assert root < domain.comm.size, f"Cannot gather on process {root} of {domain.comm.size}" + +pieces = domain.comm.gather(function_grid, root=root) + +plotter = None +if pieces is not None: + plotter = pyvista.Plotter() + plotter.open_gif("deformation.gif", fps=3) + + +def write_frame(displacements: list[np.ndarray], magnitudes: list[np.ndarray]): + """Update the values of every piece, and draw the merged grid into a single frame.""" + assert plotter is not None and pieces is not None + for piece, displacement, mag in zip(pieces, displacements, magnitudes): + piece["u"][:, : displacement.shape[1]] = displacement + piece["mag"][:] = mag + warped = pyvista.merge(pieces).warp_by_vector("u", factor=1) + warped.set_active_scalars("mag") + plotter.clear() + plotter.add_mesh(warped, show_edges=True, lighting=False, clim=[0, 10]) + plotter.update_scalar_bar_range([0, 10]) + plotter.write_frame() # - # Finally, we solve the problem over several time steps, updating the z-component of the traction @@ -214,14 +235,16 @@ def right(x): assert converged > 0, f"Solver did not converge with reason {converged}." print(f"Time step {n}, Number of iterations {num_its}, Load {T.value}") - function_grid["u"][:, : len(u)] = u.x.array.reshape(geometry.shape[0], len(u)) magnitude.interpolate(us) - warped.set_active_scalars("mag") - warped_n = function_grid.warp_by_vector(factor=1) - warped.points[:, :] = warped_n.points - warped.point_data["mag"][:] = magnitude.x.array - plotter.update_scalar_bar_range([0, 10]) - plotter.write_frame() -plotter.close() + + # The gathers are collective, so every process takes part in them, while only the process + # that received the pieces draws the frame + local_u = u.x.array.reshape(geometry.shape[0], len(u)) + displacements = domain.comm.gather(local_u.copy(), root=root) + magnitudes = domain.comm.gather(magnitude.x.array.copy(), root=root) + if displacements is not None and magnitudes is not None: + write_frame(displacements, magnitudes) +if plotter is not None: + plotter.close() # gif diff --git a/chapter2/linearelasticity_code.py b/chapter2/linearelasticity_code.py index f04b573c..f36b6393 100644 --- a/chapter2/linearelasticity_code.py +++ b/chapter2/linearelasticity_code.py @@ -143,22 +143,33 @@ def sigma(u): # We start by using Pyvista. # In previous tutorials, we have considered scalar values, while the following section considers vectors. +# Each process builds the grid of the cells it owns, the pieces are gathered on one process and +# merged into a single grid, which is then drawn, see +# [Plotting in parallel](../chapter1/fundamentals_code). + # + -# Create plotter and pyvista grid -p = pyvista.Plotter() +# Create pyvista grid topology, cell_types, geometry = plot.vtk_mesh(V) grid = pyvista.UnstructuredGrid(topology, cell_types, geometry) # Attach vector values to grid and warp grid by vector grid["u"] = uh.x.array.reshape((geometry.shape[0], 3)) -actor_0 = p.add_mesh(grid, style="wireframe", color="k") -warped = grid.warp_by_vector("u", factor=1.5) -actor_1 = p.add_mesh(warped, show_edges=True) -p.show_axes() -if not pyvista.OFF_SCREEN: - p.show() -else: - figure_as_array = p.screenshot("deflection.png") + +root = 0 +assert root < domain.comm.size, f"Cannot gather on process {root} of {domain.comm.size}" + +pieces = domain.comm.gather(grid, root=root) + +if pieces is not None: + merged_grid = pyvista.merge(pieces) + p = pyvista.Plotter() + actor_0 = p.add_mesh(merged_grid, style="wireframe", color="k") + actor_1 = p.add_mesh(merged_grid.warp_by_vector("u", factor=1.5), show_edges=True) + p.show_axes() + if not pyvista.OFF_SCREEN: + p.show() + else: + figure_as_array = p.screenshot("deflection.png") # - # We could also use Paraview for visualizing this. @@ -198,12 +209,27 @@ def sigma(u): # The first thing we notice is that we now set values for each cell, # which has a one to one correspondence with the degrees of freedom in the function space. -warped.cell_data["VonMises"] = stresses.x.petsc_vec.array -warped.set_active_scalars("VonMises") -p = pyvista.Plotter() -p.add_mesh(warped) -p.show_axes() -if not pyvista.OFF_SCREEN: - p.show() -else: - stress_figure = p.screenshot("stresses.png") +# The stresses are attached to the grid of each process before the pieces are gathered, so that +# each piece carries its own values. `petsc_vec.array` holds the cells the process owns, which are +# exactly the cells of its grid. + +# + +grid.cell_data["VonMises"] = stresses.x.petsc_vec.array +grid.set_active_scalars("VonMises") + +stress_clim = [ + domain.comm.reduce(stresses.x.petsc_vec.array.min(), op=MPI.MIN, root=root), + domain.comm.reduce(stresses.x.petsc_vec.array.max(), op=MPI.MAX, root=root), +] +stress_pieces = domain.comm.gather(grid, root=root) + +if stress_pieces is not None: + merged_stress_grid = pyvista.merge(stress_pieces) + p = pyvista.Plotter() + p.add_mesh(merged_stress_grid.warp_by_vector("u", factor=1.5), clim=stress_clim) + p.show_axes() + if not pyvista.OFF_SCREEN: + p.show() + else: + stress_figure = p.screenshot("stresses.png") +# - diff --git a/chapter2/ns_code1.py b/chapter2/ns_code1.py index 7ca88041..b6a9f319 100644 --- a/chapter2/ns_code1.py +++ b/chapter2/ns_code1.py @@ -456,6 +456,12 @@ def u_exact(x): # We have already looked at how to plot higher order functions and vector functions. # In this section we will look at how to visualize vector functions with glyphs, instead of warping the mesh. +# Each process builds the grid of the cells it owns, the pieces are gathered on one process and +# merged into a single grid, which is then drawn, see +# [Plotting in parallel](../chapter1/fundamentals_code). +# As the grids are merged, {py:meth}`glyph` places a single arrow at +# a point that is shared by two processes. + # + topology, cell_types, geometry = vtk_mesh(V) values = np.zeros((geometry.shape[0], 3), dtype=np.float64) @@ -464,23 +470,31 @@ def u_exact(x): # Create a point cloud of glyphs function_grid = pyvista.UnstructuredGrid(topology, cell_types, geometry) function_grid["u"] = values -glyphs = function_grid.glyph(orient="u", factor=0.2) # Create a pyvista-grid for the mesh tdim = mesh.topology.dim mesh.topology.create_connectivity(tdim, tdim) grid = pyvista.UnstructuredGrid(*vtk_mesh(mesh, tdim)) +root = 0 +assert root < mesh.comm.size, f"Cannot gather on process {root} of {mesh.comm.size}" + +u_pieces = mesh.comm.gather(function_grid, root=root) +mesh_pieces = mesh.comm.gather(grid, root=root) + # Create plotter -plotter = pyvista.Plotter() -plotter.add_mesh(grid, style="wireframe", color="k") -plotter.add_mesh(glyphs) -plotter.view_xy() - -if not pyvista.OFF_SCREEN: - plotter.show() -else: - fig_as_array = plotter.screenshot("glyphs.png") +if u_pieces is not None and mesh_pieces is not None: + merged_u_grid = pyvista.merge(u_pieces) + merged_grid = pyvista.merge(mesh_pieces) + plotter = pyvista.Plotter() + plotter.add_mesh(merged_grid, style="wireframe", color="k") + plotter.add_mesh(merged_u_grid.glyph(orient="u", factor=0.2)) + plotter.view_xy() + + if not pyvista.OFF_SCREEN: + plotter.show() + else: + fig_as_array = plotter.screenshot("glyphs.png") # - # ## References diff --git a/chapter3/component_bc.py b/chapter3/component_bc.py index 0d2d0847..48e92fd8 100644 --- a/chapter3/component_bc.py +++ b/chapter3/component_bc.py @@ -174,9 +174,12 @@ def sigma(u): # ## Visualization +# Each process builds the grid of the cells it owns, the pieces are gathered on one process and +# merged into a single grid, which is then drawn, see +# [Plotting in parallel](../chapter1/fundamentals_code). + # + -# Create plotter and pyvista grid -p = pyvista.Plotter() +# Create pyvista grid topology, cell_types, x = vtk_mesh(V) grid = pyvista.UnstructuredGrid(topology, cell_types, x) @@ -185,11 +188,19 @@ def sigma(u): vals = np.zeros((x.shape[0], 3)) vals[:, : len(uh)] = uh.x.array.reshape((x.shape[0], len(uh))) grid["u"] = vals -actor_0 = p.add_mesh(grid, style="wireframe", color="k") -warped = grid.warp_by_vector("u", factor=1.5) -actor_1 = p.add_mesh(warped, opacity=0.8) -p.view_xy() -if not pyvista.OFF_SCREEN: - p.show() -else: - fig_array = p.screenshot("component.png") + +root = 0 +assert root < mesh.comm.size, f"Cannot gather on process {root} of {mesh.comm.size}" + +pieces = mesh.comm.gather(grid, root=root) + +if pieces is not None: + merged_grid = pyvista.merge(pieces) + p = pyvista.Plotter() + actor_0 = p.add_mesh(merged_grid, style="wireframe", color="k") + actor_1 = p.add_mesh(merged_grid.warp_by_vector("u", factor=1.5), opacity=0.8) + p.view_xy() + if not pyvista.OFF_SCREEN: + p.show() + else: + fig_array = p.screenshot("component.png") diff --git a/chapter3/em.py b/chapter3/em.py index 0fdf59f2..a1f85aa7 100644 --- a/chapter3/em.py +++ b/chapter3/em.py @@ -226,19 +226,36 @@ # We can also visualize the subdommains using pyvista -plotter = pyvista.Plotter() +# Each process builds the grid of the cells it owns, the pieces are gathered on one process and +# merged into a single grid, which is then drawn, see +# [Plotting in parallel](../chapter1/fundamentals_code). + tdim = mesh.topology.dim mesh.topology.create_connectivity(tdim, tdim) grid = pyvista.UnstructuredGrid(*vtk_mesh(mesh, tdim)) num_local_cells = mesh.topology.index_map(tdim).size_local -grid.cell_data["Marker"] = ct.values[ct.indices < num_local_cells] +markers = ct.values[ct.indices < num_local_cells] +grid.cell_data["Marker"] = markers grid.set_active_scalars("Marker") -actor = plotter.add_mesh(grid, show_edges=True) -plotter.view_xy() -if not pyvista.OFF_SCREEN: - plotter.show() -else: - cell_tag_fig = plotter.screenshot("cell_tags.png") + +root = 0 +assert root < mesh.comm.size, f"Cannot gather on process {root} of {mesh.comm.size}" + +marker_clim = [ + mesh.comm.reduce(markers.min(), op=MPI.MIN, root=root), + mesh.comm.reduce(markers.max(), op=MPI.MAX, root=root), +] +marker_pieces = mesh.comm.gather(grid, root=root) + +if marker_pieces is not None: + marker_grid = pyvista.merge(marker_pieces) + plotter = pyvista.Plotter() + actor = plotter.add_mesh(marker_grid, show_edges=True, clim=marker_clim) + plotter.view_xy() + if not pyvista.OFF_SCREEN: + plotter.show() + else: + cell_tag_fig = plotter.screenshot("cell_tags.png") # Next, we define the discontinous functions for the permeability $\mu$ and current $J_z$ using the `MeshTags` as in [Defining material parameters through subdomains](./subdomains) @@ -300,17 +317,26 @@ # We now plot the magnetic potential $A_z$ and the magnetic field $B$. We start by creating a new plotter # + -plotter = pyvista.Plotter() - Az_grid = pyvista.UnstructuredGrid(*vtk_mesh(V)) Az_grid.point_data["A_z"] = A_z.x.array Az_grid.set_active_scalars("A_z") -warp = Az_grid.warp_by_scalar("A_z", factor=1e7) -actor = plotter.add_mesh(warp, show_edges=True) -if not pyvista.OFF_SCREEN: - plotter.show() -else: - Az_fig = plotter.screenshot("Az.png") + +Az_clim = [ + mesh.comm.reduce(A_z.x.array.min(), op=MPI.MIN, root=root), + mesh.comm.reduce(A_z.x.array.max(), op=MPI.MAX, root=root), +] +Az_pieces = mesh.comm.gather(Az_grid, root=root) + +if Az_pieces is not None: + merged_Az_grid = pyvista.merge(Az_pieces) + plotter = pyvista.Plotter() + actor = plotter.add_mesh( + merged_Az_grid.warp_by_scalar("A_z", factor=1e7), show_edges=True, clim=Az_clim + ) + if not pyvista.OFF_SCREEN: + plotter.show() + else: + Az_fig = plotter.screenshot("Az.png") # - # ## Visualizing the magnetic field @@ -320,28 +346,35 @@ # We connect the vector field with the midpoint by using `pyvista.PolyData`. # + -plotter = pyvista.Plotter() -plotter.set_position([0, 0, 5]) - -# We include ghosts cells as we access all degrees of freedom (including ghosts) on each process +# We use the cells owned by the process, so that no arrow is drawn twice where two processes +# meet, and gather the clouds along with the grid of the mesh top_imap = mesh.topology.index_map(mesh.topology.dim) -num_cells = top_imap.size_local + top_imap.num_ghosts +num_cells = top_imap.size_local mesh.topology.create_connectivity(mesh.topology.dim, mesh.topology.dim) midpoints = compute_midpoints( mesh, mesh.topology.dim, np.arange(num_cells, dtype=np.int32) ) -num_dofs = W.dofmap.index_map.size_local + W.dofmap.index_map.num_ghosts +num_dofs = W.dofmap.index_map.size_local assert num_cells == num_dofs +bs = W.dofmap.index_map_bs values = np.zeros((num_dofs, 3), dtype=np.float64) -values[:, : mesh.geometry.dim] = B.x.array.real.reshape(num_dofs, W.dofmap.index_map_bs) +values[:, : mesh.geometry.dim] = B.x.array.real[: num_dofs * bs].reshape(num_dofs, bs) cloud = pyvista.PolyData(midpoints) cloud["B"] = values -glyphs = cloud.glyph("B", factor=2e6) -actor = plotter.add_mesh(grid, style="wireframe", color="k") -actor2 = plotter.add_mesh(glyphs) - -if not pyvista.OFF_SCREEN: - plotter.show() -else: - B_fig = plotter.screenshot("B.png") + +cloud_pieces = mesh.comm.gather(cloud, root=root) +mesh_pieces = mesh.comm.gather(grid, root=root) + +if cloud_pieces is not None and mesh_pieces is not None: + merged_cloud = pyvista.merge(cloud_pieces) + merged_grid = pyvista.merge(mesh_pieces) + plotter = pyvista.Plotter() + plotter.set_position([0, 0, 5]) + actor = plotter.add_mesh(merged_grid, style="wireframe", color="k") + actor2 = plotter.add_mesh(merged_cloud.glyph("B", factor=2e6)) + + if not pyvista.OFF_SCREEN: + plotter.show() + else: + B_fig = plotter.screenshot("B.png") diff --git a/chapter3/multiple_dirichlet.py b/chapter3/multiple_dirichlet.py index 9edf23d8..1e1d6f4f 100644 --- a/chapter3/multiple_dirichlet.py +++ b/chapter3/multiple_dirichlet.py @@ -121,16 +121,31 @@ def u_exact(x): # To visualize the solution, run the script with in a Jupyter notebook with `off_screen=False` or as a python script with `off_screen=True`. # + +# Each process builds the grid of the cells it owns, the pieces are gathered on one process and +# merged into a single grid, which is then drawn, see +# [Plotting in parallel](../chapter1/fundamentals_code). + pyvista_cells, cell_types, geometry = vtk_mesh(V) grid = pyvista.UnstructuredGrid(pyvista_cells, cell_types, geometry) grid.point_data["u"] = uh.x.array grid.set_active_scalars("u") -plotter = pyvista.Plotter() -plotter.add_text("uh", position="upper_edge", font_size=14, color="black") -plotter.add_mesh(grid, show_edges=True) -plotter.view_xy() -if not pyvista.OFF_SCREEN: - plotter.show() -else: - figure = plotter.screenshot("multiple_dirichlet.png") +root = 0 +assert root < mesh.comm.size, f"Cannot gather on process {root} of {mesh.comm.size}" + +clim = [ + mesh.comm.reduce(uh.x.array.min(), op=MPI.MIN, root=root), + mesh.comm.reduce(uh.x.array.max(), op=MPI.MAX, root=root), +] +pieces = mesh.comm.gather(grid, root=root) + +if pieces is not None: + merged_grid = pyvista.merge(pieces) + plotter = pyvista.Plotter() + plotter.add_text("uh", position="upper_edge", font_size=14, color="black") + plotter.add_mesh(merged_grid, show_edges=True, clim=clim) + plotter.view_xy() + if not pyvista.OFF_SCREEN: + plotter.show() + else: + figure = plotter.screenshot("multiple_dirichlet.png") diff --git a/chapter3/neumann_dirichlet_code.py b/chapter3/neumann_dirichlet_code.py index 397d0202..0a4ce56d 100644 --- a/chapter3/neumann_dirichlet_code.py +++ b/chapter3/neumann_dirichlet_code.py @@ -186,17 +186,32 @@ def boundary_D(x): # To look at the actual solution, run the script as a python script with `off_screen=True` or as a Jupyter notebook with `off_screen=False` # + +# Each process builds the grid of the cells it owns, the pieces are gathered on one process and +# merged into a single grid, which is then drawn, see +# [Plotting in parallel](../chapter1/fundamentals_code). + pyvista_cells, cell_types, geometry = vtk_mesh(V) grid = pyvista.UnstructuredGrid(pyvista_cells, cell_types, geometry) grid.point_data["u"] = uh.x.array grid.set_active_scalars("u") -plotter = pyvista.Plotter() -plotter.add_text("uh", position="upper_edge", font_size=14, color="black") -plotter.add_mesh(grid, show_edges=True) -plotter.view_xy() - -if not pyvista.OFF_SCREEN: - plotter.show() -else: - figure = plotter.screenshot("neumann_dirichlet.png") +root = 0 +assert root < mesh.comm.size, f"Cannot gather on process {root} of {mesh.comm.size}" + +clim = [ + mesh.comm.reduce(uh.x.array.min(), op=MPI.MIN, root=root), + mesh.comm.reduce(uh.x.array.max(), op=MPI.MAX, root=root), +] +pieces = mesh.comm.gather(grid, root=root) + +if pieces is not None: + merged_grid = pyvista.merge(pieces) + plotter = pyvista.Plotter() + plotter.add_text("uh", position="upper_edge", font_size=14, color="black") + plotter.add_mesh(merged_grid, show_edges=True, clim=clim) + plotter.view_xy() + + if not pyvista.OFF_SCREEN: + plotter.show() + else: + figure = plotter.screenshot("neumann_dirichlet.png") diff --git a/chapter3/robin_neumann_dirichlet.py b/chapter3/robin_neumann_dirichlet.py index a068e7c6..029ddacf 100644 --- a/chapter3/robin_neumann_dirichlet.py +++ b/chapter3/robin_neumann_dirichlet.py @@ -267,19 +267,34 @@ def type(self): uh = problem.solve() # Visualize solution +# Each process builds the grid of the cells it owns, the pieces are gathered on one process and +# merged into a single grid, which is then drawn, see +# [Plotting in parallel](../chapter1/fundamentals_code). + pyvista_cells, cell_types, geometry = vtk_mesh(V) grid = pyvista.UnstructuredGrid(pyvista_cells, cell_types, geometry) grid.point_data["u"] = uh.x.array grid.set_active_scalars("u") -plotter = pyvista.Plotter() -plotter.add_text("uh", position="upper_edge", font_size=14, color="black") -plotter.add_mesh(grid, show_edges=True) -plotter.view_xy() -if not pyvista.OFF_SCREEN: - plotter.show() -else: - figure = plotter.screenshot("robin_neumann_dirichlet.png") +root = 0 +assert root < mesh.comm.size, f"Cannot gather on process {root} of {mesh.comm.size}" + +clim = [ + mesh.comm.reduce(uh.x.array.min(), op=MPI.MIN, root=root), + mesh.comm.reduce(uh.x.array.max(), op=MPI.MAX, root=root), +] +pieces = mesh.comm.gather(grid, root=root) + +if pieces is not None: + merged_grid = pyvista.merge(pieces) + plotter = pyvista.Plotter() + plotter.add_text("uh", position="upper_edge", font_size=14, color="black") + plotter.add_mesh(merged_grid, show_edges=True, clim=clim) + plotter.view_xy() + if not pyvista.OFF_SCREEN: + plotter.show() + else: + figure = plotter.screenshot("robin_neumann_dirichlet.png") # - # ## Verification diff --git a/chapter3/subdomains.py b/chapter3/subdomains.py index 2815ec04..1b6e1bc4 100644 --- a/chapter3/subdomains.py +++ b/chapter3/subdomains.py @@ -127,6 +127,10 @@ def Omega_1(x): ) uh = problem.solve() +# Each process builds the grid of the cells it owns, the pieces are gathered on one process and +# merged into a single grid, which is then drawn, see +# [Plotting in parallel](../chapter1/fundamentals_code). + # Filter out ghosted cells tdim = mesh.topology.dim num_cells_local = mesh.topology.index_map(tdim).size_local @@ -140,25 +144,46 @@ def Omega_1(x): mesh, tdim, np.arange(num_cells_local, dtype=np.int32) ) -p = pyvista.Plotter(window_size=[800, 800]) grid = pyvista.UnstructuredGrid(topology, cell_types, x) grid.cell_data["Marker"] = marker grid.set_active_scalars("Marker") -p.add_mesh(grid, show_edges=True) -if pyvista.OFF_SCREEN: - figure = p.screenshot("subdomains_structured.png") -p.show() + +root = 0 +assert root < mesh.comm.size, f"Cannot gather on process {root} of {mesh.comm.size}" + +marker_clim = [ + mesh.comm.reduce(marker.min(), op=MPI.MIN, root=root), + mesh.comm.reduce(marker.max(), op=MPI.MAX, root=root), +] +marker_pieces = mesh.comm.gather(grid, root=root) + +if marker_pieces is not None: + marker_grid = pyvista.merge(marker_pieces) + p = pyvista.Plotter(window_size=[800, 800]) + p.add_mesh(marker_grid, show_edges=True, clim=marker_clim) + if pyvista.OFF_SCREEN: + figure = p.screenshot("subdomains_structured.png") + p.show() # - -p2 = pyvista.Plotter(window_size=[800, 800]) grid_uh = pyvista.UnstructuredGrid(*vtk_mesh(V)) grid_uh.point_data["u"] = uh.x.array.real grid_uh.set_active_scalars("u") -p2.add_mesh(grid_uh, show_edges=True) -if not pyvista.OFF_SCREEN: - p2.show() -else: - figure = p2.screenshot("subdomains_structured2.png") + +u_clim = [ + mesh.comm.reduce(uh.x.array.real.min(), op=MPI.MIN, root=root), + mesh.comm.reduce(uh.x.array.real.max(), op=MPI.MAX, root=root), +] +u_pieces = mesh.comm.gather(grid_uh, root=root) + +if u_pieces is not None: + merged_u_grid = pyvista.merge(u_pieces) + p2 = pyvista.Plotter(window_size=[800, 800]) + p2.add_mesh(merged_u_grid, show_edges=True, clim=u_clim) + if not pyvista.OFF_SCREEN: + p2.show() + else: + figure = p2.screenshot("subdomains_structured2.png") # We clearly observe different behavior in the two regions, which both have the same Dirichlet boundary condition on the left side, where $x=0$. @@ -341,22 +366,40 @@ def create_mesh(mesh, cell_type, prune_z=False): topology, cell_types, x = vtk_mesh(mesh, tdim) grid = pyvista.UnstructuredGrid(topology, cell_types, x) num_local_cells = mesh.topology.index_map(tdim).size_local -grid.cell_data["Marker"] = ct.values[ct.indices < num_local_cells] +markers = ct.values[ct.indices < num_local_cells] +grid.cell_data["Marker"] = markers grid.set_active_scalars("Marker") -p = pyvista.Plotter(window_size=[800, 800]) -p.add_mesh(grid, show_edges=True) -if not pyvista.OFF_SCREEN: - p.show() -else: - figure = p.screenshot("subdomains_unstructured.png") +marker_clim = [ + mesh.comm.reduce(markers.min(), op=MPI.MIN, root=root), + mesh.comm.reduce(markers.max(), op=MPI.MAX, root=root), +] +marker_pieces = mesh.comm.gather(grid, root=root) + +if marker_pieces is not None: + marker_grid = pyvista.merge(marker_pieces) + p = pyvista.Plotter(window_size=[800, 800]) + p.add_mesh(marker_grid, show_edges=True, clim=marker_clim) + if not pyvista.OFF_SCREEN: + p.show() + else: + figure = p.screenshot("subdomains_unstructured.png") # - grid_uh = pyvista.UnstructuredGrid(*vtk_mesh(V)) grid_uh.point_data["u"] = uh.x.array.real grid_uh.set_active_scalars("u") -p2 = pyvista.Plotter(window_size=[800, 800]) -p2.add_mesh(grid_uh, show_edges=True) -if not pyvista.OFF_SCREEN: - p2.show() -else: - p2.screenshot("unstructured_u.png") + +u_clim = [ + mesh.comm.reduce(uh.x.array.real.min(), op=MPI.MIN, root=root), + mesh.comm.reduce(uh.x.array.real.max(), op=MPI.MAX, root=root), +] +u_pieces = mesh.comm.gather(grid_uh, root=root) + +if u_pieces is not None: + merged_u_grid = pyvista.merge(u_pieces) + p2 = pyvista.Plotter(window_size=[800, 800]) + p2.add_mesh(merged_u_grid, show_edges=True, clim=u_clim) + if not pyvista.OFF_SCREEN: + p2.show() + else: + p2.screenshot("unstructured_u.png") diff --git a/chapter4/newton-solver.py b/chapter4/newton-solver.py index 97c9e184..b9781ebc 100644 --- a/chapter4/newton-solver.py +++ b/chapter4/newton-solver.py @@ -333,12 +333,28 @@ def u_exact(x): error_max = domain.comm.allreduce(np.max(np.abs(uh.x.array - u_D.x.array)), op=MPI.MAX) PETSc.Sys.Print(f"Error_max: {error_max:.2e}") +# Each process builds the grid of the cells it owns, the pieces are gathered on one process and +# merged into a single grid, which is then drawn, see +# [Plotting in parallel](../chapter1/fundamentals_code). + u_topology, u_cell_types, u_geometry = dolfinx.plot.vtk_mesh(V) u_grid = pyvista.UnstructuredGrid(u_topology, u_cell_types, u_geometry) u_grid.point_data["u"] = uh.x.array.real u_grid.set_active_scalars("u") -u_plotter = pyvista.Plotter() -u_plotter.add_mesh(u_grid, show_edges=True) -u_plotter.view_xy() -if not pyvista.OFF_SCREEN: - u_plotter.show() + +root = 0 +assert root < domain.comm.size, f"Cannot gather on process {root} of {domain.comm.size}" + +u_clim = [ + domain.comm.reduce(uh.x.array.real.min(), op=MPI.MIN, root=root), + domain.comm.reduce(uh.x.array.real.max(), op=MPI.MAX, root=root), +] +pieces = domain.comm.gather(u_grid, root=root) + +if pieces is not None: + merged_u_grid = pyvista.merge(pieces) + u_plotter = pyvista.Plotter() + u_plotter.add_mesh(merged_u_grid, show_edges=True, clim=u_clim) + u_plotter.view_xy() + if not pyvista.OFF_SCREEN: + u_plotter.show() From 19019f4dfbba96117d78c0e4e9ff07173ef33e8a Mon Sep 17 00:00:00 2001 From: jorgensd Date: Wed, 16 Sep 2026 22:55:11 +0200 Subject: [PATCH 2/5] Add extra note --- chapter1/fundamentals_code.py | 16 ++++++++++++---- 1 file changed, 12 insertions(+), 4 deletions(-) diff --git a/chapter1/fundamentals_code.py b/chapter1/fundamentals_code.py index 018b26f5..934e93b1 100644 --- a/chapter1/fundamentals_code.py +++ b/chapter1/fundamentals_code.py @@ -327,10 +327,18 @@ # single grid with {py:func}`pyvista.merge`. # {py:func}`dolfinx.plot.vtk_mesh` uses only the cells that the process *owns*, so a cell that is # shared between two processes is not drawn twice. -# `merge` welds the points that the pieces have in common back together, so the merged grid behaves -# like one built on a single process: filters that follow the field from one cell into the next, -# such as {py:meth}`streamlines`, do not stop at the partition -# boundaries. +# `merge` welds the points that the pieces have in common back together, so we get a single mesh +# that is drawn with a single color bar, and that filters can be applied to as a whole. +# +# ```{admonition} Points on a partition boundary +# :class: tip +# The processes on either side of a partition boundary compute the coordinates of the points they +# share independently, so these can differ in the last bits, and `merge` then keeps both copies. +# The picture is the same either way, but a filter that follows the field from one cell into the +# next, such as {py:meth}`streamlines`, would stop at such a +# point. Merge with a tolerance to weld them as well: +# `pieces[0].merge(pieces[1:], tolerance=1e-14)`. +# ``` # # The process that collects the pieces is chosen with `root`, which has to be one of the processes # we run on. From 3c1d4fdaa975c534ef3f6f84b0d9b8f8986b6a60 Mon Sep 17 00:00:00 2001 From: jorgensd Date: Thu, 17 Sep 2026 05:56:10 +0200 Subject: [PATCH 3/5] Formatting --- chapter1/fundamentals_code.py | 5 ++++- chapter2/hyperelasticity.py | 2 +- 2 files changed, 5 insertions(+), 2 deletions(-) diff --git a/chapter1/fundamentals_code.py b/chapter1/fundamentals_code.py index 934e93b1..3214465d 100644 --- a/chapter1/fundamentals_code.py +++ b/chapter1/fundamentals_code.py @@ -421,7 +421,10 @@ if u_pieces is not None: plotter2 = pyvista.Plotter() plotter2.add_mesh( - merged_u_grid.warp_by_scalar(), show_edges=True, show_scalar_bar=True, clim=u_clim + merged_u_grid.warp_by_scalar(), + show_edges=True, + show_scalar_bar=True, + clim=u_clim, ) if not pyvista.OFF_SCREEN: plotter2.show() diff --git a/chapter2/hyperelasticity.py b/chapter2/hyperelasticity.py index 2559a438..c00fc90e 100644 --- a/chapter2/hyperelasticity.py +++ b/chapter2/hyperelasticity.py @@ -210,7 +210,7 @@ def right(x): def write_frame(displacements: list[np.ndarray], magnitudes: list[np.ndarray]): - """Update the values of every piece, and draw the merged grid into a single frame.""" + """Update the values of each piece, and draw the merged grid as a single frame.""" assert plotter is not None and pieces is not None for piece, displacement, mag in zip(pieces, displacements, magnitudes): piece["u"][:, : displacement.shape[1]] = displacement From 5a895b45f4dcd11643168933567a39a40dd09eeb Mon Sep 17 00:00:00 2001 From: jorgensd Date: Thu, 17 Sep 2026 05:57:01 +0200 Subject: [PATCH 4/5] Update AMR --- chapter2/amr.py | 5 +---- 1 file changed, 1 insertion(+), 4 deletions(-) diff --git a/chapter2/amr.py b/chapter2/amr.py index 0ad3c892..45f795e2 100644 --- a/chapter2/amr.py +++ b/chapter2/amr.py @@ -291,10 +291,7 @@ def write_frame(plotter: pyvista.Plotter | None, uh_r: dolfinx.fem.Function): if plotter is None or pieces is None or curved_pieces is None: return plotter.clear() - plotter.add_mesh( - pyvista.merge(curved_pieces), - show_edges=False, - ) + plotter.add_mesh(pyvista.merge(curved_pieces), show_edges=False) plotter.add_mesh(pyvista.merge(pieces), style="wireframe", color="black") plotter.view_xy() plotter.write_frame() From 9f09bfafcdb76eeb15775a7f7a596cc9c86a5077 Mon Sep 17 00:00:00 2001 From: jorgensd Date: Thu, 17 Sep 2026 06:00:20 +0200 Subject: [PATCH 5/5] Add pyc to gitignore --- .gitignore | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/.gitignore b/.gitignore index 91fd7f7c..94e1a3b0 100644 --- a/.gitignore +++ b/.gitignore @@ -9,4 +9,5 @@ _build *.png *.pvtu *.msh -*.bp \ No newline at end of file +*.bp +*.pyc \ No newline at end of file