mesh_smoother_test.C

Go to the documentation of this file.

#include <libmesh/boundary_info.h>
#include <libmesh/elem.h>
#include <libmesh/enum_elem_type.h>
#include <libmesh/function_base.h>
#include <libmesh/libmesh.h>
#include <libmesh/mesh_generation.h>
#include <libmesh/mesh_modification.h>
#include <libmesh/mesh_smoother_laplace.h>
#include <libmesh/mesh_smoother_vsmoother.h>
#include <libmesh/mesh_tools.h>
#include <libmesh/node.h>
#include <libmesh/reference_elem.h>
#include <libmesh/replicated_mesh.h>
#include <libmesh/system.h> // LIBMESH_HAVE_SOLVER define
#include "libmesh/face_tri.h"
#include "libmesh/cell_prism.h"
#include "libmesh/utility.h"
#include "libmesh/enum_to_string.h"
#include "libmesh/parallel_ghost_sync.h"

#include "test_comm.h"
#include "libmesh_cppunit.h"

#include <random>

namespace
{
using namespace libMesh;
using Utility::pow;

// Distortion function that doesn't distort boundary nodes
// 2D only, use for LaplaceMeshSmoother
class DistortSquare : public FunctionBase<Real>
{
  std::unique_ptr<FunctionBase<Real>> clone() const override
  {
    return std::make_unique<DistortSquare>();
  }

  Real operator()(const Point &, const Real = 0.) override
  {
    libmesh_not_implemented();
  } // scalar-only API

  // Skew inward based on a cubic function
  void operator()(const Point & p, const Real, DenseVector<Real> & output)
  {
    output.resize(3);
    const Real eta = 2 * p(0) - 1;
    const Real zeta = 2 * p(1) - 1;
    output(0) = p(0) + (pow<3>(eta) - eta) * p(1) * (1 - p(1));
    output(1) = p(1) + (pow<3>(zeta) - zeta) * p(0) * (1 - p(0));
    output(2) = 0;
  }
};

// Distortion function supporting 1D, 2D, and 3D
// Use for VariationalMeshSmoother
class DistortHyperCube : public FunctionBase<Real>
{
public:
  DistortHyperCube(const unsigned int dim) : _dim(dim) {}

private:
  std::unique_ptr<FunctionBase<Real>> clone() const override
  {
    return std::make_unique<DistortHyperCube>(_dim);
  }

  Real operator()(const Point &, const Real = 0.) override { libmesh_not_implemented(); }

  void operator()(const Point & p, const Real, DenseVector<Real> & output) override
  {
    output.resize(3);
    output.zero();

    // Count how many coordinates are exactly on the boundary (0 or 1)
    std::array<bool, 3> is_on_boundary = {false, false, false};
    for (unsigned int i = 0; i < _dim; ++i)
      if (std::abs(p(i)) < TOLERANCE || std::abs(p(i) - 1.) < TOLERANCE)
        is_on_boundary[i] = true;

    // Distort only those directions not fixed on the boundary
    for (unsigned int i = 0; i < _dim; ++i)
      {
        if (!is_on_boundary[i]) // only distort free dimensions
          {
            // This value controls the strength of the distortion
            Real modulation = 0.3;
            Real xi = 2. * p(i) - 1.;
            for (unsigned int j = 0; j < _dim; ++j)
              {
                if (j != i)
                  {
                    Real pj = std::clamp(p(j), 0., 1.);         // ensure numeric safety
                    modulation *= (pj - 0.5) * (pj - 0.5) * 4.; // quadratic bump centered at 0.5
                  }
              }
            const auto delta = (pow<3>(xi) - xi) * modulation;
            // Check for delta = 0 to make sure we perturb every point
            output(i) = (delta > TOLERANCE) ? p(i) + delta : 1.03 * p(i);
          }
        else
          output(i) = p(i); // dimension on boundary remains unchanged
      }
  }

  const unsigned int _dim;
};

class SquareToParallelogram : public FunctionBase<Real>
{
  std::unique_ptr<FunctionBase<Real>> clone() const override
  {
    return std::make_unique<SquareToParallelogram>();
  }

  Real operator()(const Point &, const Real = 0.) override
  {
    libmesh_not_implemented();
  } // scalar-only API

  // Has the effect that a square, meshed into right triangles with diagonals
  // rising from lower-left to upper-right, is transformed into a left-leaning
  // parallelogram of eqilateral triangles.
  void operator()(const Point & p, const Real, DenseVector<Real> & output)
  {
    output.resize(3);
    output(0) = p(0) - 0.5 * p(1);
    output(1) = 0.5 * std::sqrt(Real(3)) * p(1);
    output(2) = p(2);
  }
};

class ParallelogramToSquare : public FunctionBase<Real>
{
  std::unique_ptr<FunctionBase<Real>> clone() const override
  {
    return std::make_unique<ParallelogramToSquare>();
  }

  Real operator()(const Point &, const Real = 0.) override
  {
    libmesh_not_implemented();
  } // scalar-only API

  // Has the effect that a left-leaning parallelogram of equilateral triangles is transformed
  // into a square of right triangle with diagonals rising from lower-left to upper-right.
  // This is the inversion of the SquareToParallelogram mapping.
  void operator()(const Point & p, const Real, DenseVector<Real> & output)
  {
    output.resize(3);
    output(0) = p(0) + p(1) / std::sqrt(Real(3));
    output(1) = (2. / std::sqrt(Real(3))) * p(1);
    output(2) = p(2);
  }
};

class CubeToParallelepiped : public FunctionBase<Real>
{
  std::unique_ptr<FunctionBase<Real>> clone() const override
  {
    return std::make_unique<CubeToParallelepiped>();
  }

  Real operator()(const Point &, const Real = 0.) override
  {
    libmesh_not_implemented();
  } // scalar-only API

  // Has the effect that a cube, meshed into right triangular prisms with diagonals
  // rising from lower-right to upper-left, is transformed into a right-leaning
  // parallelepiped of equilateral triangular prisms. Additionally, the z direction
  // is scaled to ensure element height to base aspect ratios match the target element.
  // Without the correct aspect ratio, the smoothed mesh nodes are off and the
  // distortion_is assertions fail.
  void operator()(const Point & p, const Real, DenseVector<Real> & output)
  {
    output.resize(3);
    output(0) = p(0) + 0.5 * p(1);
    output(1) = 0.5 * std::sqrt(3) * p(1);
    // Adjusting z to get the triangular area to element height aspect ratio To
    // match the target element
    output(2) = p(2) * 0.25 * std::sqrt(3.);
  }
};

class ParallelepipedToCube : public FunctionBase<Real>
{
  std::unique_ptr<FunctionBase<Real>> clone() const override
  {
    return std::make_unique<ParallelepipedToCube>();
  }

  Real operator()(const Point &, const Real = 0.) override
  {
    libmesh_not_implemented();
  } // scalar-only API

  // Has the effect that a right-leaning parallelepiped of equilaterial triangular
  // prisms with is transformed into a cube of right triangular prisms with
  // diagonals rising from lower-right to upper-left. Additionally, the z direction
  // is scaled to ensure element heights align with the values expected by the
  // distortion_is function.
  // This is the inversion of the CubeToParallelepiped mapping.
  void operator()(const Point & p, const Real, DenseVector<Real> & output)
  {
    output.resize(3);
    output(0) = p(0) - p(1) / std::sqrt(3.);
    output(1) = (2. / std::sqrt(3)) * p(1);
    // Undoing the aspect ratio adjustment to get the z divisions to work with
    // distortion_is
    output(2) = p(2) * 4. / std::sqrt(3.);
  }
};
}

using namespace libMesh;

class MeshSmootherTest : public CppUnit::TestCase
{
public:
  LIBMESH_CPPUNIT_TEST_SUITE(MeshSmootherTest);

#if LIBMESH_DIM > 2
  CPPUNIT_TEST(testLaplaceQuad);
  CPPUNIT_TEST(testLaplaceTri);
#if defined(LIBMESH_ENABLE_VSMOOTHER) && defined(LIBMESH_HAVE_SOLVER)
  CPPUNIT_TEST(testVariationalEdge2);
  CPPUNIT_TEST(testVariationalEdge3);
  CPPUNIT_TEST(testVariationalEdge3MultipleSubdomains);

  CPPUNIT_TEST(testVariationalQuad4);
  CPPUNIT_TEST(testVariationalQuad4MultipleSubdomains);
  CPPUNIT_TEST(testVariationalQuad8);
  CPPUNIT_TEST(testVariationalQuad9);
  CPPUNIT_TEST(testVariationalQuad9MultipleSubdomains);
  CPPUNIT_TEST(testVariationalQuad4Tangled);

  CPPUNIT_TEST(testVariationalTri3);
  CPPUNIT_TEST(testVariationalTri6);
  CPPUNIT_TEST(testVariationalTri6MultipleSubdomains);

  CPPUNIT_TEST(testVariationalHex8);
  CPPUNIT_TEST(testVariationalHex20);
  CPPUNIT_TEST(testVariationalHex27);
  CPPUNIT_TEST(testVariationalHex27Tangled);
  CPPUNIT_TEST(testVariationalHex27MultipleSubdomains);

  CPPUNIT_TEST(testVariationalPyramid5);
  CPPUNIT_TEST(testVariationalPyramid13);
  CPPUNIT_TEST(testVariationalPyramid14);
  CPPUNIT_TEST(testVariationalPyramid18);
  CPPUNIT_TEST(testVariationalPyramid18MultipleSubdomains);

  CPPUNIT_TEST(testVariationalPrism6);
  CPPUNIT_TEST(testVariationalPrism15);
  CPPUNIT_TEST(testVariationalPrism18);
  CPPUNIT_TEST(testVariationalPrism20);
  CPPUNIT_TEST(testVariationalPrism21);
  CPPUNIT_TEST(testVariationalPrism21MultipleSubdomains);

#if defined(LIBMESH_HAVE_GZSTREAM)
  CPPUNIT_TEST(testVariationalMixed2D);
  CPPUNIT_TEST(testVariationalMixed3D);
#endif

#endif // LIBMESH_ENABLE_VSMOOTHER
#endif

  CPPUNIT_TEST_SUITE_END();

public:
  void setUp() {}

  void tearDown() {}

  void testLaplaceSmoother(ReplicatedMesh & mesh, MeshSmoother & smoother, ElemType type)
  {
    LOG_UNIT_TEST;

    constexpr unsigned int n_elems_per_side = 4;

    MeshTools::Generation::build_square(
        mesh, n_elems_per_side, n_elems_per_side, 0., 1., 0., 1., type);

    // Move it around so we have something that needs smoothing
    DistortSquare ds;
    MeshTools::Modification::redistribute(mesh, ds);

    // Assert the distortion is as expected
    auto center_distortion_is =
        [](const Node & node, int d, bool distortion, Real distortion_tol = TOLERANCE) {
          const Real r = node(d);
          const Real R = r * n_elems_per_side;
          CPPUNIT_ASSERT_GREATER(-TOLERANCE * TOLERANCE, r);
          CPPUNIT_ASSERT_GREATER(-TOLERANCE * TOLERANCE, 1 - r);

          // If we're at the boundaries we should *never* be distorted
          if (std::abs(node(0)) < TOLERANCE * TOLERANCE ||
              std::abs(node(0) - 1) < TOLERANCE * TOLERANCE)
            {
              const Real R1 = node(1) * n_elems_per_side;
              CPPUNIT_ASSERT_LESS(TOLERANCE * TOLERANCE, std::abs(R1 - std::round(R1)));
              return true;
            }

          if (std::abs(node(1)) < TOLERANCE * TOLERANCE ||
              std::abs(node(1) - 1) < TOLERANCE * TOLERANCE)
            {
              const Real R0 = node(0) * n_elems_per_side;
              CPPUNIT_ASSERT_LESS(TOLERANCE * TOLERANCE, std::abs(R0 - std::round(R0)));

              return true;
            }

          // If we're at the center we're fine
          if (std::abs(r - 0.5) < TOLERANCE * TOLERANCE)
            return true;

          return ((std::abs(R - std::round(R)) > distortion_tol) == distortion);
        };

    for (auto node : mesh.node_ptr_range())
      {
        CPPUNIT_ASSERT(center_distortion_is(*node, 0, true));
        CPPUNIT_ASSERT(center_distortion_is(*node, 1, true));
      }

    // Enough iterations to mostly fix us up.  Laplace seems to be at 1e-3
    // tolerance by iteration 6, so hopefully everything is there on any
    // system by 8.
    for (unsigned int i = 0; i != 8; ++i)
      smoother.smooth();

    // Make sure we're not too distorted anymore.
    for (auto node : mesh.node_ptr_range())
      {
        CPPUNIT_ASSERT(center_distortion_is(*node, 0, false, 1e-3));
        CPPUNIT_ASSERT(center_distortion_is(*node, 1, false, 1e-3));
      }
  }

  void testVariationalSmoother(ReplicatedMesh & mesh,
                               VariationalMeshSmoother & smoother,
                               const ElemType type,
                               const bool multiple_subdomains = false,
                               const bool tangle_mesh=false,
                               const Real tangle_damping_factor=1.0)
  {
    LOG_UNIT_TEST;

    if (multiple_subdomains && tangle_mesh)
      libmesh_not_implemented_msg(
          "Arbitrary mesh tangling with multiple subdomains is not supported.");

    // Get mesh dimension, determine whether type is triangular
    const auto * ref_elem = &(ReferenceElem::get(type));
    const auto dim = ref_elem->dim();
    const bool type_is_tri = Utility::enum_to_string(type).compare(0, 3, "TRI") == 0;
    const bool type_is_pyramid = Utility::enum_to_string(type).compare(0, 7, "PYRAMID") == 0;
    const bool type_is_prism = Utility::enum_to_string(type).compare(0, 5, "PRISM") == 0;

    // Used fewer elems for higher order types, as extra midpoint nodes will add
    // enough complexity
    unsigned int n_elems_per_side = 4 / Elem::type_to_default_order_map[type];

    switch (dim)
      {
        case 1:
          MeshTools::Generation::build_line(mesh, n_elems_per_side, 0., 1., type);
          break;
        case 2:
          MeshTools::Generation::build_square(
              mesh, n_elems_per_side, n_elems_per_side, 0., 1., 0., 1., type);
          break;

        case 3:
          MeshTools::Generation::build_cube(mesh,
                                            n_elems_per_side,
                                            n_elems_per_side,
                                            n_elems_per_side,
                                            0.,
                                            1.,
                                            0.,
                                            1.,
                                            0.,
                                            1.,
                                            type);
          break;

        default:
          libmesh_error_msg("Unsupported dimension " << dim);
      }

    // Move it around so we have something that needs smoothing
    DistortHyperCube dh(dim);
    MeshTools::Modification::redistribute(mesh, dh);

    const auto & boundary_info = mesh.get_boundary_info();

    if (tangle_mesh)
      {
        // We tangle the mesh by (partially) swapping the locations of 2 nodes.
        // If the n-dimentional hypercube is represented as a grid with
        // (undistorted) divisions occuring at integer numbers, we will
        // (partially) swap the nodes nearest p1 = (1,1,1) and p2 = (2,1,2).

        // Define p1 and p2 given the integer indices
        // Undistorted elem side length
        const Real dr = 1. / (n_elems_per_side * Elem::type_to_default_order_map[type]);
        const Point p1 = Point(dr, dim > 1 ? dr : 0, dim > 2 ? dr : 0);
        const Point p2 = Point(2. * dr, dim > 1 ? dr : 0, dim > 2 ? 2. * dr : 0);

        // ids and distances of mesh nodes nearest p1 and p2
        dof_id_type node1_id = DofObject::invalid_id;
        Real dist1_closest = std::numeric_limits<Real>::max();
        dof_id_type node2_id = DofObject::invalid_id;
        Real dist2_closest = std::numeric_limits<Real>::max();

        for (const auto * node : mesh.local_node_ptr_range())
          {
            const auto dist1 = (*node - p1).norm();
            if (dist1 < dist1_closest)
              {
                dist1_closest = dist1;
                node1_id = node->id();
              }

            const auto dist2 = (*node - p2).norm();
            if (dist2 < dist2_closest)
              {
                dist2_closest = dist2;
                node2_id = node->id();
              }
          }

        // parallel communication
        unsigned int node1_rank;
        mesh.comm().minloc(dist1_closest, node1_rank);
        mesh.comm().broadcast(node1_id, node1_rank);

        unsigned int node2_rank;
        mesh.comm().minloc(dist2_closest, node2_rank);
        mesh.comm().broadcast(node2_id, node2_rank);

        // modify the nodes
        if ((node1_id != DofObject::invalid_id) && (node2_id != DofObject::invalid_id))
          {
            libmesh_assert(node1_id != node2_id);

            auto & node1 = mesh.node_ref(node1_id);
            auto & node2 = mesh.node_ref(node2_id);

            const Point diff = node1 - node2;
            // Note that a damping factor of 1 will swap the nodes and one
            // of 0.5 will move the nodes to the mean of their original
            // locations.
            node1 -= tangle_damping_factor * diff;
            node2 += tangle_damping_factor * diff;
          }

        SyncNodalPositions sync_object(mesh);
        Parallel::sync_dofobject_data_by_id (mesh.comm(), mesh.nodes_begin(), mesh.nodes_end(), sync_object);

        // Make sure the mesh is tangled
        smoother.setup(); // Need this to create the system we are about to access
        const auto & unsmoothed_info = smoother.get_mesh_info();
        CPPUNIT_ASSERT(unsmoothed_info.mesh_is_tangled);
      }

    // Add multiple subdomains if requested
    std::unordered_map<subdomain_id_type, Real> distorted_subdomain_volumes;
    if (multiple_subdomains)
      {
        // Modify the subdomain ids in an interesting way
        for (auto * elem : mesh.active_element_ptr_range())
          {
            unsigned int subdomain_id = 0;
            for (const auto d : make_range(dim))
              if (elem->vertex_average()(d) > 0.5)
                ++subdomain_id;
            elem->subdomain_id() += subdomain_id;
          }

        // This loop should NOT be combined with the one above because we need to
        // finish checking and updating subdomain ids for all elements before
        // recording the final subdomain volumes.
        for (auto * elem : mesh.active_element_ptr_range())
          distorted_subdomain_volumes[elem->subdomain_id()] += elem->volume();
      }

    // Get the mesh order
    const auto & elem_orders = mesh.elem_default_orders();
    libmesh_error_msg_if(elem_orders.size() != 1,
                         "The variational smoother cannot be used for mixed-order meshes!");
    // For pyramids, the factor 2 accounts for the account that cubes of side
    // length s are divided into pyramids of base side length s and height s/2.
    // The factor 4 lets us catch triangular face midpoint nodes for PYRAMID18 elements.
    const auto scale_factor = *elem_orders.begin() * (type_is_pyramid ? 2 * 4 : 1);

    // Function to assert the node distortion is as expected
    auto node_distortion_is = [&n_elems_per_side, &dim, &boundary_info, &scale_factor, &type_is_prism](
                             const Node & node, bool distortion, Real distortion_tol = TOLERANCE) {
      // Get boundary ids associated with the node
      std::vector<boundary_id_type> boundary_ids;
      boundary_info.boundary_ids(&node, boundary_ids);

      // This tells us what type of node we are: internal, sliding, or fixed
      const auto num_dofs = dim - boundary_ids.size();
      /*
       * The following cases of num_dofs are possible, ASSUMING all boundaries
       * are non-overlapping
       * 3D: 3-0,     3-1,     3-2,     3-3
       *   = 3        2        1        0
       *     internal sliding, sliding, fixed
       * 2D: 2-0,     2-1,     2-2
       *   = 2        1        0
       *     internal sliding, fixed
       * 1D: 1-0,     1-1
       *   = 1        0
       *     internal fixed
       *
       * We expect that R is an integer in [0, n_elems_per_side] for
       * num_dofs of the node's cooridinantes, while the remaining coordinates
       * are fixed to the boundary with value 0 or 1. In other words, at LEAST
       * dim - num_dofs coordinantes should be 0 or 1.
       */

      std::size_t num_zero_or_one = 0;

      bool distorted = false;
      for (const auto d : make_range(dim))
        {
          const Real r = node(d);
          const Real R = r * n_elems_per_side * scale_factor;
          CPPUNIT_ASSERT_GREATER(-distortion_tol * distortion_tol, r);
          CPPUNIT_ASSERT_GREATER(-distortion_tol * distortion_tol, 1 - r);

          bool d_distorted = std::abs(R - std::round(R)) > distortion_tol;
          if (type_is_prism && (scale_factor == 3))
            {
              // Adjustment required for triangular face nodes of PRISM20/21 elements.
              // These elements have fe_order 3. This takes care of nodes occuring at
              // thirds, but not halves.
              const Real R_prism = R / scale_factor * 2;
              const bool d_distorted_prism =
                  std::abs(R_prism - std::round(R_prism)) > distortion_tol;
              d_distorted &= d_distorted_prism;
            }
          distorted |= d_distorted;
          num_zero_or_one += (absolute_fuzzy_equals(r, 0.) || absolute_fuzzy_equals(r, 1.));
        }

      CPPUNIT_ASSERT_GREATEREQUAL(dim - num_dofs, num_zero_or_one);

      // We can never expect a fixed node to be distorted
      if (num_dofs == 0)
        // if (num_dofs < dim)
        return true;
      return distorted == distortion;
    };

    // Make sure our DistortHyperCube transformation has distorted the mesh
    for (auto node : mesh.node_ptr_range())
      CPPUNIT_ASSERT(node_distortion_is(*node, true));

    // Transform the square mesh of triangles to a parallelogram mesh of
    // triangles. This will allow the Variational Smoother to smooth the mesh
    // to the optimal case of equilateral triangles
    if (type_is_tri)
      {
        // Transform the square mesh of triangles to a parallelogram mesh of
        // triangles. This will allow the Variational Smoother to smooth the mesh
        // to the optimal case of equilateral triangles
        SquareToParallelogram stp;
        MeshTools::Modification::redistribute(mesh, stp);
      }
    else if (type_is_prism)
      {
        // Transform the cube mesh of prisms to a parallelepiped mesh of
        // prisms. This will allow the Variational Smoother to smooth the mesh
        // to the optimal case of equilateral triangular prisms
        CubeToParallelepiped ctp;
        MeshTools::Modification::redistribute(mesh, ctp);
      }

    smoother.smooth();

    if (type_is_tri)
      {
        // Transform the parallelogram mesh back to a square mesh. In the case of
        // the Variational Smoother, equilateral triangules will be transformed
        // into right triangules that align with the original undistorted mesh.
        ParallelogramToSquare pts;
        MeshTools::Modification::redistribute(mesh, pts);
      }
    else if (type_is_prism)
      {
        // Transform the parallelepiped mesh back to a cube mesh. In the case of
        // the Variational Smoother, equilateral triangular prisms will be transformed
        // into right triangular prisms that align with the original undistorted mesh.
        ParallelepipedToCube ptc;
        MeshTools::Modification::redistribute(mesh, ptc);
      }

    if (multiple_subdomains)
      {
        // Make sure the subdomain volumes didn't change. Although this doesn't
        // guarantee that the subdomain didn't change, it is a good indicator
        // and is friedly to sliding subdomain boundary nodes.
        std::unordered_map<subdomain_id_type, Real> smoothed_subdomain_volumes;
        for (auto * elem : mesh.active_element_ptr_range())
          smoothed_subdomain_volumes[elem->subdomain_id()] += elem->volume();

        for (const auto & pair : distorted_subdomain_volumes)
          {
            const auto & subdomain_id = pair.first;
            CPPUNIT_ASSERT(
                relative_fuzzy_equals(libmesh_map_find(distorted_subdomain_volumes, subdomain_id),
                                      libmesh_map_find(smoothed_subdomain_volumes, subdomain_id),
                                      TOLERANCE));
          }
      }
    else
      {
        // Make sure we're not too distorted anymore
        std::set<dof_id_type> nodes_checked;
        for (const auto * elem : mesh.active_element_ptr_range())
          {
            for (const auto local_node_id : make_range(elem->n_nodes()))
              {
                const auto & node = elem->node_ref(local_node_id);
                if (nodes_checked.find(node.id()) != nodes_checked.end())
                  continue;

                nodes_checked.insert(node.id());

                // Check for special case where pyramidal base-to-apex midpoint
                // nodes are not smoothed to the actual midpoints.
                if (type_is_pyramid)
                  {
                    if (local_node_id > 8 && local_node_id < 13)
                      {
                        // Base-Apex midpoint nodes
                        //
                        // Due to the nature of the pyramidal target element and the
                        // cubic nature of the test mesh, for a dilation weight o
                        // 0.5, smoothed base-apex midpoint nodes are smoothed to
                        // the point base + x * (apex - base) instead of
                        // base + 0.5 (apex - base), where x is in (0,1) and
                        // depends on the number of nodes in the element.
                        // We hard-code this check below.

                        const auto & base = elem->node_ref(local_node_id - 9);
                        const auto & apex = elem->node_ref(4);
                        const Real x = (type == PYRAMID18) ? 0.569332 : 0.549876;

                        CPPUNIT_ASSERT(node.relative_fuzzy_equals(base + x * (apex - base), 1e-3));
                        continue;
                      }
                    else if (local_node_id > 13)
                      {
                        // Triangular face midpoint nodes
                        //
                        // Due to the nature of the pyramidal target element and the
                        // cubic nature of the test mesh, for a dilation weight o
                        // 0.5, smoothed triangular face midpoint nodes are
                        // smoothed a weighted average of the constituent
                        // vertices instead of an unweighted average.
                        // We hard-code this check below.
                        const auto & base1 = elem->node_ref(local_node_id - 14);
                        const auto & base2 = elem->node_ref((local_node_id - 13) % 4);
                        const auto & apex = elem->node_ref(4);

                        const auto node_approx = (0.3141064847 * base1 +
                                                  0.3141064847 * base2 +
                                                  0.3717870306 * apex);
                        CPPUNIT_ASSERT(node.relative_fuzzy_equals(node_approx, 1e-3));
                        continue;
                      }
                  }

                CPPUNIT_ASSERT(node_distortion_is(node, false, 1e-2));

              }
          }
      }
  }

  // Function that distorts and smooths a mesh and checks that the result is
  // equal to the original "gold" mesh
  void testVariationalSmootherRegression(const ReplicatedMesh & gold_mesh)
  {
    // Make a copy of the gold mesh to distort and then smooth
    ReplicatedMesh mesh(gold_mesh);

    // Move it around so we have something that needs smoothing
    DistortHyperCube dh(mesh.mesh_dimension());
    MeshTools::Modification::redistribute(mesh, dh);

    // Function to assert the distortion is as expected
    auto mesh_distortion_is = [](const ReplicatedMesh & gold_mesh,
                                 const ReplicatedMesh & mesh,
                                 const bool distortion,
                                 Real distortion_tol = TOLERANCE) {

      CPPUNIT_ASSERT(gold_mesh.n_nodes() ==  mesh.n_nodes());
      for (const auto node_id : make_range(gold_mesh.n_nodes()))
        {
          const auto & gold_node = gold_mesh.node_ref(node_id);
          const auto & node = mesh.node_ref(node_id);
          for (const auto d : make_range(gold_mesh.mesh_dimension()))
            {
              const bool d_distorted = std::abs(gold_node(d) - node(d)) > distortion_tol;
              if (d_distorted)
                // mesh is distorted from gold_mesh
                return distortion;
            }
        }

      // mesh is not distorted from gold_mesh
      return !distortion;
    };

    // Make sure the mesh has been distorted
    CPPUNIT_ASSERT(mesh_distortion_is(gold_mesh, mesh, true));

    // Turn off subdomain boundary preservation because DistortHyperCube
    // does not preserve subdomain boundaries
    // Also set dilation coefficient to 0 because the reference volume of the
    // distorted mesh won't be equal to the reference volume of the smoothed
    // mesh and will smooth to a slightly different solution.
    VariationalMeshSmoother smoother(mesh, 0.0, false);
    smoother.smooth();

    // Make sure the mesh has been smoothed to the gold mesh
    CPPUNIT_ASSERT(mesh_distortion_is(gold_mesh, mesh, false));
  }

  void testLaplaceQuad()
  {
    ReplicatedMesh mesh(*TestCommWorld);
    LaplaceMeshSmoother laplace(mesh);

    testLaplaceSmoother(mesh, laplace, QUAD4);
  }

  void testLaplaceTri()
  {
    ReplicatedMesh mesh(*TestCommWorld);
    LaplaceMeshSmoother laplace(mesh);

    testLaplaceSmoother(mesh, laplace, TRI3);
  }

#ifdef LIBMESH_ENABLE_VSMOOTHER
  void testVariationalEdge2()
  {
    ReplicatedMesh mesh(*TestCommWorld);
    VariationalMeshSmoother variational(mesh);

    testVariationalSmoother(mesh, variational, EDGE2);
  }

  void testVariationalEdge3()
  {
    ReplicatedMesh mesh(*TestCommWorld);
    VariationalMeshSmoother variational(mesh);

    testVariationalSmoother(mesh, variational, EDGE3);
  }

  void testVariationalEdge3MultipleSubdomains()
  {
    ReplicatedMesh mesh(*TestCommWorld);
    VariationalMeshSmoother variational(mesh);

    testVariationalSmoother(mesh, variational, EDGE3, true);
  }

  void testVariationalQuad4()
  {
    ReplicatedMesh mesh(*TestCommWorld);
    VariationalMeshSmoother variational(mesh);

    testVariationalSmoother(mesh, variational, QUAD4);
  }

  void testVariationalQuad4MultipleSubdomains()
  {
    ReplicatedMesh mesh(*TestCommWorld);
    VariationalMeshSmoother variational(mesh);

    testVariationalSmoother(mesh, variational, QUAD4, true);
  }

  void testVariationalQuad8()
  {
    ReplicatedMesh mesh(*TestCommWorld);
    VariationalMeshSmoother variational(mesh);

    testVariationalSmoother(mesh, variational, QUAD8);
  }

  void testVariationalQuad9()
  {
    ReplicatedMesh mesh(*TestCommWorld);
    VariationalMeshSmoother variational(mesh);

    testVariationalSmoother(mesh, variational, QUAD9);
  }

  void testVariationalQuad9MultipleSubdomains()
  {
    ReplicatedMesh mesh(*TestCommWorld);
    VariationalMeshSmoother variational(mesh);

    testVariationalSmoother(mesh, variational, QUAD9, true);
  }

  void testVariationalQuad4Tangled()
  {
    ReplicatedMesh mesh(*TestCommWorld);
    VariationalMeshSmoother variational(mesh);

    testVariationalSmoother(mesh, variational, QUAD4, false, true, 0.65);
  }

  void testVariationalTri3()
  {
    ReplicatedMesh mesh(*TestCommWorld);
    VariationalMeshSmoother variational(mesh);

    testVariationalSmoother(mesh, variational, TRI3);
  }

  void testVariationalTri6()
  {
    ReplicatedMesh mesh(*TestCommWorld);
    VariationalMeshSmoother variational(mesh);

    testVariationalSmoother(mesh, variational, TRI6);
  }

  void testVariationalTri6MultipleSubdomains()
  {
    ReplicatedMesh mesh(*TestCommWorld);
    VariationalMeshSmoother variational(mesh);

    testVariationalSmoother(mesh, variational, TRI3, true);
  }

  void testVariationalHex8()
  {
    ReplicatedMesh mesh(*TestCommWorld);
    VariationalMeshSmoother variational(mesh);

    testVariationalSmoother(mesh, variational, HEX8);
  }

  void testVariationalHex20()
  {
    ReplicatedMesh mesh(*TestCommWorld);
    VariationalMeshSmoother variational(mesh);

    testVariationalSmoother(mesh, variational, HEX20);
  }

  void testVariationalHex27()
  {
    ReplicatedMesh mesh(*TestCommWorld);
    VariationalMeshSmoother variational(mesh);

    testVariationalSmoother(mesh, variational, HEX27);
  }

  void testVariationalHex27MultipleSubdomains()
  {
    ReplicatedMesh mesh(*TestCommWorld);
    VariationalMeshSmoother variational(mesh);

    testVariationalSmoother(mesh, variational, HEX27, true);
  }

  void testVariationalHex27Tangled()
  {
    ReplicatedMesh mesh(*TestCommWorld);
    VariationalMeshSmoother variational(mesh);

    testVariationalSmoother(mesh, variational, HEX27, false, true, 0.55);
  }

  void testVariationalPyramid5()
  {
    ReplicatedMesh mesh(*TestCommWorld);
    VariationalMeshSmoother variational(mesh);

    testVariationalSmoother(mesh, variational, PYRAMID5);
  }

  void testVariationalPyramid13()
  {
    ReplicatedMesh mesh(*TestCommWorld);
    VariationalMeshSmoother variational(mesh);

    testVariationalSmoother(mesh, variational, PYRAMID13);
  }

  void testVariationalPyramid14()
  {
    ReplicatedMesh mesh(*TestCommWorld);
    VariationalMeshSmoother variational(mesh);

    testVariationalSmoother(mesh, variational, PYRAMID14);
  }

  void testVariationalPyramid18()
  {
    ReplicatedMesh mesh(*TestCommWorld);
    VariationalMeshSmoother variational(mesh);

    testVariationalSmoother(mesh, variational, PYRAMID18);
  }

  void testVariationalPyramid18MultipleSubdomains()
  {
    ReplicatedMesh mesh(*TestCommWorld);
    VariationalMeshSmoother variational(mesh);

    testVariationalSmoother(mesh, variational, PYRAMID18, true);
  }

  void testVariationalPrism6()
  {
    ReplicatedMesh mesh(*TestCommWorld);
    VariationalMeshSmoother variational(mesh);

    testVariationalSmoother(mesh, variational, PRISM6);
  }

  void testVariationalPrism15()
  {
    ReplicatedMesh mesh(*TestCommWorld);
    VariationalMeshSmoother variational(mesh);

    testVariationalSmoother(mesh, variational, PRISM15);
  }

  void testVariationalPrism18()
  {
    ReplicatedMesh mesh(*TestCommWorld);
    VariationalMeshSmoother variational(
        mesh, 0.5, true, TOLERANCE * TOLERANCE, 10 * TOLERANCE * TOLERANCE);

    testVariationalSmoother(mesh, variational, PRISM18);
  }

  void testVariationalPrism20()
  {
    ReplicatedMesh mesh(*TestCommWorld);
    VariationalMeshSmoother variational(
        mesh, 0.5, true, TOLERANCE * TOLERANCE, 10 * TOLERANCE * TOLERANCE);

    testVariationalSmoother(mesh, variational, PRISM20);
  }

  void testVariationalPrism21()
  {
    ReplicatedMesh mesh(*TestCommWorld);
    VariationalMeshSmoother variational(
        mesh, 0.5, true, TOLERANCE * TOLERANCE, 10 * TOLERANCE * TOLERANCE);

    testVariationalSmoother(mesh, variational, PRISM21);
  }

  void testVariationalPrism21MultipleSubdomains()
  {
    ReplicatedMesh mesh(*TestCommWorld);
    mesh.allow_renumbering(false);
    VariationalMeshSmoother variational(mesh);

    testVariationalSmoother(mesh, variational, PRISM21, true);
  }

  void testVariationalMixed2D()
  {
    ReplicatedMesh mesh(*TestCommWorld);
    mesh.read("meshes/quad4_tri3_smoothed.xda.gz");

    testVariationalSmootherRegression(mesh);
  }

  void testVariationalMixed3D()
  {
    ReplicatedMesh mesh(*TestCommWorld);
    mesh.read("meshes/hex8_prism6_smoothed.xda.gz");

    testVariationalSmootherRegression(mesh);
  }
#endif // LIBMESH_ENABLE_VSMOOTHER
};

CPPUNIT_TEST_SUITE_REGISTRATION(MeshSmootherTest);