constraint_manager.h

Go to the documentation of this file.

// SPDX-FileCopyrightText: © 2025 PRISMS Center at the University of Michigan
// SPDX-License-Identifier: GNU Lesser General Public Version 2.1
 
#pragma once
 
#include <deal.II/base/mg_level_object.h>
#include <deal.II/dofs/dof_handler.h>
#include <deal.II/fe/mapping.h>
#include <deal.II/lac/affine_constraints.h>
#include <deal.II/multigrid/mg_constrained_dofs.h>
 
#include <prismspf/core/dof_manager.h>
#include <prismspf/core/field_attributes.h>
#include <prismspf/core/solve_block.h>
#include <prismspf/core/type_enums.h>
#include <prismspf/core/types.h>
 
#include <prismspf/user_inputs/constraint_parameters.h>
#include <prismspf/user_inputs/spatial_discretization.h>
 
#include <prismspf/config.h>
 
PRISMS_PF_BEGIN_NAMESPACE
 
// Forward declaration to avoid circular dependency.
template <unsigned int dim, unsigned int degree, typename number>
class PDEOperatorBase;
 
// TODO (fractalsbyx): The following snippet is from dealii. Using this (by
// pre-constructing the needed maps and functions) may be cleaner than how dirichlet are
// currently done.
/* template <int dim, int spacedim, typename number>
   void
  interpolate_boundary_values(
    const Mapping<dim, spacedim>    &mapping,
    const DoFHandler<dim, spacedim> &dof,
    const std::map<types::boundary_id, const Function<spacedim, number> *>
                              &function_map,
    AffineConstraints<number> &constraints,
    const ComponentMask       &component_mask = {}); */
 
// TODO (fractalsbyx): I don't think we actually need 'change' constraints for newton
// solves at all. Verify this and remove them.

template <unsigned int dim, unsigned int degree, typename number>
class ConstraintManager
{
public:
  ConstraintManager() = default;

  void
  init(const BoundaryParameters<dim>              &_boundary_parameters,
       const SpatialDiscretization<dim>           &_spatial_discretization,
       const DoFManager<dim, degree>              &_dof_manager,
       const PDEOperatorBase<dim, degree, number> &_pde_operator);

  [[nodiscard]] std::vector<const dealii::AffineConstraints<number> *>
  get_field_constraints() const;

  [[nodiscard]] const dealii::AffineConstraints<number> &
  get_field_constraint(Types::Index index) const;

  [[nodiscard]] const std::array<dealii::AffineConstraints<number>, 2> &
  get_generic_constraints() const;

  [[nodiscard]] const dealii::AffineConstraints<number> &
  get_generic_constraint(unsigned int rank) const;
 
  //-------------------------------------------------------------------------------------

  [[nodiscard]] const std::vector<std::vector<dealii::AffineConstraints<number>>> &
  get_mg_field_constraints_levels() const;

  [[nodiscard]] std::vector<const dealii::AffineConstraints<number> *>
  get_mg_field_constraints(unsigned int relative_level) const;

  [[nodiscard]] const dealii::AffineConstraints<number> &
  get_mg_field_constraint(Types::Index index, unsigned int relative_level) const;

  [[nodiscard]] const std::vector<std::array<dealii::AffineConstraints<number>, 2>> &
  get_mg_generic_constraints_levels() const;

  [[nodiscard]] const std::array<dealii::AffineConstraints<number>, 2> &
  get_mg_generic_constraints(unsigned int relative_level) const;

  [[nodiscard]] const dealii::AffineConstraints<number> &
  get_mg_generic_constraint(unsigned int rank, unsigned int relative_level) const;

  void
  reinit(const std::vector<FieldAttributes> &field_attributes);

  void
  update_time_dependent_constraints(const std::vector<FieldAttributes> &field_attributes);

  static const std::array<dealii::ComponentMask, dim> vector_component_mask;
  static const dealii::ComponentMask                  scalar_empty_mask;
 
private:
  void
  make_constraints_for_single_field(dealii::AffineConstraints<number> &constraint,
                                    const dealii::DoFHandler<dim>     &dof_handler,
                                    const FieldConstraints<dim>       &_field_constraints,
                                    TensorRank                         tensor_rank,
                                    Types::Index                       field_index,
                                    unsigned int relative_level = -1);

  void
  make_bc_constraints(dealii::AffineConstraints<number> &constraint,
                      const dealii::DoFHandler<dim>     &dof_handler,
                      const FieldConstraints<dim>       &boundary_condition,
                      TensorRank                         tensor_rank,
                      Types::Index                       field_index);

  void
  make_one_boundary_constraint(dealii::AffineConstraints<number> &_constraints,
                               unsigned int                       boundary_id,
                               unsigned int                       component,
                               Condition                          boundary_type,
                               const dealii::DoFHandler<dim>     &dof_handler,
                               TensorRank                         tensor_rank,
                               Types::Index                       field_index) const;

  void
  make_natural_constraints() const;

  void
  make_dirichlet_constraints(dealii::AffineConstraints<number> &_constraints,
                             const dealii::DoFHandler<dim>     &dof_handler,
                             const unsigned int                &boundary_id,
                             const unsigned int                &field_index,
                             const bool                        &is_vector_field,
                             const dealii::ComponentMask       &mask) const;

  void
  make_periodic_constraints(dealii::AffineConstraints<number> &_constraints,
                            const dealii::DoFHandler<dim>     &dof_handler,
                            const unsigned int                &boundary_id,
                            const dealii::ComponentMask       &mask) const;

  const BoundaryParameters<dim> *boundary_parameters = nullptr;

  const SpatialDiscretization<dim> *spatial_discretization = nullptr;

  const DoFManager<dim, degree> *dof_manager = nullptr;

  const PDEOperatorBase<dim, degree, number> *pde_operator = nullptr;

  std::vector<dealii::AffineConstraints<number>> field_constraints;

  std::array<dealii::AffineConstraints<number>, 2> generic_constraints;

  std::vector<std::vector<dealii::AffineConstraints<number>>> mg_field_constraints;

  std::vector<std::array<dealii::AffineConstraints<number>, 2>> mg_generic_constraints;
};
 
template <unsigned int dim, unsigned int degree, typename number>
inline const std::array<dealii::ComponentMask, dim>
  ConstraintManager<dim, degree, number>::vector_component_mask = []()
{
  std::array<dealii::ComponentMask, dim> masks {};
  for (unsigned int i = 0; i < dim; ++i)
    {
      dealii::ComponentMask temp_mask(dim, false);
      temp_mask.set(i, true);
      masks.at(i) = temp_mask;
    }
  return masks;
}();
 
template <unsigned int dim, unsigned int degree, typename number>
inline const dealii::ComponentMask
  ConstraintManager<dim, degree, number>::scalar_empty_mask {};
 
PRISMS_PF_END_NAMESPACE