variational_smoother_constraint.h

Go to the documentation of this file.

// The libMesh Finite Element Library.
// Copyright (C) 2002-2025 Benjamin S. Kirk, John W. Peterson, Roy H. Stogner

// This library is free software; you can redistribute it and/or
// modify it under the terms of the GNU Lesser General Public
// License as published by the Free Software Foundation; either
// version 2.1 of the License, or (at your option) any later version.

// This library is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
// Lesser General Public License for more details.

// You should have received a copy of the GNU Lesser General Public
// License along with this library; if not, write to the Free Software
// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA

#ifndef LIBMESH_VARIATIONAL_SMOOTHER_CONSTRAINT_H
#define LIBMESH_VARIATIONAL_SMOOTHER_CONSTRAINT_H

// Local Includes
#include "libmesh/system.h"
#include "libmesh/dof_map.h"

// C++ includes
#include <variant>

namespace libMesh
{

// Forward declarations
class PointConstraint;
class LineConstraint;
class PlaneConstraint;
class InvalidConstraint;

using ConstraintVariant = std::variant<PointConstraint, LineConstraint,
                                       PlaneConstraint, InvalidConstraint>;

class PointConstraint
{

public:
  PointConstraint() = default;

  PointConstraint(const Point &point, const Real &tol = TOLERANCE * TOLERANCE);

  bool operator<(const PointConstraint &other) const;

  bool operator==(const PointConstraint &other) const;

  bool contains_point(const PointConstraint &p) const { return *this == p; }

  ConstraintVariant intersect(const ConstraintVariant &other) const;

  const Point &point() const { return _point; }

  const Real &tol() const { return _tol; }

private:
  // Life is easier if we don't make this const
  Point _point;

  Real _tol;
};

class LineConstraint
{
public:
  LineConstraint() = default;

  LineConstraint(const Point &point, const Point &direction,
                 const Real &tol = TOLERANCE * TOLERANCE);

  bool operator<(const LineConstraint &other) const;

  bool operator==(const LineConstraint &other) const;

  bool contains_point(const PointConstraint &p) const;

  bool is_parallel(const LineConstraint &l) const;

  bool is_parallel(const PlaneConstraint &p) const;

  ConstraintVariant intersect(const ConstraintVariant &other) const;

  const Point &point() const { return _point; }

  const Point &direction() const { return _direction; }

  const Real &tol() const { return _tol; }

private:
  // Life is easier if we don't make these const
  Point _point;

  Point _direction;

  Real _tol;
};

class PlaneConstraint
{

public:
  PlaneConstraint() = default;

  PlaneConstraint(const Point &point, const Point &normal,
                  const Real &tol = TOLERANCE * TOLERANCE);

  bool operator<(const PlaneConstraint &other) const;

  bool operator==(const PlaneConstraint &other) const;

  bool contains_point(const PointConstraint &p) const;

  bool contains_line(const LineConstraint &l) const;

  bool is_parallel(const PlaneConstraint &p) const;

  bool is_parallel(const LineConstraint &l) const;

  ConstraintVariant intersect(const ConstraintVariant &other) const;

  const Point &point() const { return _point; }

  const Point &normal() const { return _normal; }

  const Real &tol() const { return _tol; }

private:
  // Life is easier if we don't make these const
  Point _point;

  Point _normal;

  Real _tol;
};

class InvalidConstraint
{

public:
  InvalidConstraint()
    : _err_msg("We should never get here! The InvalidConstraint object should be "
               "detected and replaced with a valid ConstraintVariant prior to calling "
               "any class methods.")
  {
  }

  ConstraintVariant intersect(const ConstraintVariant &) const {
    libmesh_assert_msg(false, _err_msg);
    return *this;
  }

  bool contains_point(const PointConstraint &) const {
    libmesh_assert_msg(false, _err_msg);
    return false;
  }

private:
  std::string _err_msg;
};

inline ConstraintVariant intersect_constraints(const ConstraintVariant &a,
                                               const ConstraintVariant &b) {
  // std::visit applies the visitor v (a Callable that can be called with any
  // combination of types from Variants) to the active value inside a
  // std::Variant. This circumvents the issue that the literal ConstraintVariant
  // type does not have a method called 'intersect' (but the types defining
  // ConstraintVariant do)
  return std::visit(
      [](const auto &lhs, const auto &rhs) -> ConstraintVariant {
        return lhs.intersect(rhs);
      },
      a, b);
}

class VariationalSmootherConstraint : public System::Constraint
{
private:

  System & _sys;

  const bool _preserve_subdomain_boundaries;

  void fix_node(const Node & node);

  void constrain_node_to_plane(const Node & node, const Point & ref_normal_vec);

  void constrain_node_to_line(const Node & node, const Point & line_vec);

  static void find_nodal_or_face_neighbors(
      const MeshBase & mesh,
      const Node & node,
      const std::unordered_map<dof_id_type, std::vector<const Elem *>> & nodes_to_elem_map,
      std::vector<const Node *> & neighbors);

  static bool nodes_share_boundary_id(
      const Node & boundary_node,
      const Node & neighbor_node,
      const Elem & containing_elem,
      const BoundaryInfo & boundary_info);

  static std::set<std::set<const Node *>>
  get_neighbors_for_subdomain_constraint(
      const MeshBase &mesh, const Node &node, const subdomain_id_type sub_id,
      const std::unordered_map<dof_id_type, std::vector<const Elem *>>
          &nodes_to_elem_map);

  static std::set<std::set<const Node *>> get_neighbors_for_boundary_constraint(
      const MeshBase &mesh, const Node &node,
      const std::unordered_set<dof_id_type> &boundary_node_ids,
      const BoundaryInfo &boundary_info,
      const std::unordered_map<dof_id_type, std::vector<const Elem *>>
          &nodes_to_elem_map);

  static ConstraintVariant determine_constraint(
      const Node &node, const unsigned int dim,
      const std::set<std::set<const Node *>> &side_grouped_boundary_neighbors);

  void impose_constraint(const Node &node, const ConstraintVariant &constraint);

public:
  VariationalSmootherConstraint(System & sys, const bool & preserve_subdomain_boundaries);

  virtual ~VariationalSmootherConstraint() override;

  virtual void constrain() override;
};


} // namespace libMesh

#endif // LIBMESH_VARIATIONAL_SMOOTHER_CONSTRAINT_H