plot_modes Subroutine

subroutine plot_modes(s, X, lambda, phi, names, scale, filename)
  !! Write one PNG, SVG or PDF (chosen by the extension of filename) with a subplot per
  !! mode: undeformed (black), deformed (green), nodal displacements (blue)
  use fplotlib, only: figure, subplot, plot, plot3d, scatter, scatter3d, title, &
                    & set_aspect, axis, tight_layout, savefig
  type(shapefunc), intent(in) :: s
  real(dp), intent(in) :: X(:,:), lambda(:), phi(:,:)
  real(dp), intent(in) :: scale   !! largest displacement relative to the element size
  character(len=*), intent(in) :: names(:), filename
  real(dp) :: U(size(X,1), size(X,2)), h
  integer, allocatable :: ev(:,:), path(:)
  integer :: ndof, rows, cols, m, e, i, nv, ne
  character(len=120) :: label

  ndof = size(lambda)
  cols = ceiling(sqrt(real(ndof)))
  rows = (ndof + cols - 1)/cols
  h = maxval(maxval(X, 2) - minval(X, 2))
  ev = edge_vertices(s%elm)
  nv = maxval(ev)
  ne = max(s%order - 1, 0)   ! interior nodes per edge
  call figure(figsize=[4.5_dp*cols, 3.5_dp*rows])

  do m = 1, ndof
    call subplot(rows, cols, m)
    U = X + scale*h*reshape(phi(:,m), shape(X))/maxval(abs(phi(:,m)))
    write (label, '("$\lambda_{",i0,"}$=",es9.2,2x,a)') m, lambda(m), trim(names(m))
    call title(trim(label), fontsize=8.0_dp)

    ! 1D: undeformed at y = 0, deformed at y = h/2
    if (s%dim == 1) then
      call draw2d(reshape([X, 0*X], [2, size(X,2)], order=[2,1]), &
                & reshape([U, 0*U + h/2], [2, size(U,2)], order=[2,1]))
    else if (s%dim == 2) then
      call draw2d(X, U)
    else
      do e = 1, size(ev, 2)
        path = edge_path(e)
        call plot3d(X(1,path), X(2,path), X(3,path), "k-")
        call plot3d(U(1,path), U(2,path), U(3,path), "g-")
      end do
      do i = 1, size(X, 2)
        call plot3d([X(1,i), U(1,i)], [X(2,i), U(2,i)], [X(3,i), U(3,i)], "b-", lw=0.5_dp)
      end do
      call scatter3d(X(1,:), X(2,:), X(3,:), c="k", s=10.0_dp)
      call scatter3d(U(1,:), U(2,:), U(3,:), c="k", s=10.0_dp)
    end if
  end do
  do m = ndof + 1, rows*cols
    call subplot(rows, cols, m)
    call axis("off")
  end do
  call tight_layout()   ! trim the white border
  call savefig(filename)

contains

function edge_path(e) result(p)
  !! Nodes along edge e: first vertex, edge nodes, second vertex
  ! Allocatable result: ifx 2026.1.1 miscompiles an automatic result sized by the host's ne
  ! when draw2d calls this function (segfault, or wrong nodes with -check all).
  integer, intent(in) :: e
  integer, allocatable :: p(:)
  integer :: k

  allocate (p(ne + 2))
  p(1) = ev(1,e)
  p(ne+2) = ev(2,e)
  p(2:ne+1) = [(nv + (e - 1)*ne + k, k = 1, ne)]
end function edge_path

subroutine draw2d(A, C)
  !! Edges and nodes of the undeformed (A) and deformed (C) element in the plane
  real(dp), intent(in) :: A(:,:), C(:,:)
  integer :: e, i

  do e = 1, size(ev, 2)
    path = edge_path(e)
    call plot(A(1,path), A(2,path), "k-")
    call plot(C(1,path), C(2,path), "g-")
  end do
  do i = 1, size(A, 2)
    call plot([A(1,i), C(1,i)], [A(2,i), C(2,i)], "b-", lw=0.5_dp)
  end do
  call scatter(A(1,:), A(2,:), c="k", s=10.0_dp)
  call scatter(C(1,:), C(2,:), c="k", s=10.0_dp)
  call set_aspect(1.0_dp)
end subroutine draw2d

end subroutine plot_modes