EqualValueBoundaryConstraint.C

Go to the documentation of this file.

//* This file is part of the MOOSE framework
//* https://mooseframework.inl.gov
//*
//* All rights reserved, see COPYRIGHT for full restrictions
//* https://github.com/idaholab/moose/blob/master/COPYRIGHT
//*
//* Licensed under LGPL 2.1, please see LICENSE for details
//* https://www.gnu.org/licenses/lgpl-2.1.html
 
// MOOSE includes
#include "EqualValueBoundaryConstraint.h"
#include "DisplacedProblem.h"
#include "MooseMesh.h"
 
registerMooseObject("MooseApp", EqualValueBoundaryConstraint);
 
InputParameters
EqualValueBoundaryConstraint::validParams()
{
  InputParameters params = NodalConstraint::validParams();
  params.addClassDescription(
      "Constraint for enforcing that variables on each side of a boundary are equivalent.");
  params.addParam<unsigned int>(
      "primary",
      "The ID of the primary node. If no ID is provided, first node of secondary set is chosen.");
  params.addParam<Point>("primary_node_coord", "Coordinates of the primary node to locate.");
  params.addParam<std::vector<unsigned int>>("secondary_node_ids", "The IDs of the secondary node");
  params.addParam<BoundaryName>("secondary", "The boundary ID associated with the secondary side");
  params.addRequiredParam<Real>("penalty", "The penalty used for the boundary term");
  return params;
}
 
EqualValueBoundaryConstraint::EqualValueBoundaryConstraint(const InputParameters & parameters)
  : NodalConstraint(parameters), _penalty(getParam<Real>("penalty"))
{
  updateConstrainedNodes();
}
 
void
EqualValueBoundaryConstraint::meshChanged()
{
  updateConstrainedNodes();
}
 
void
EqualValueBoundaryConstraint::updateConstrainedNodes()
{
  populateSecondaryNodes();
  pickPrimaryNode();
  ghostPrimary();
}
 
void
EqualValueBoundaryConstraint::populateSecondaryNodes()
{
  // user provided nothing
  if (!isParamValid("secondary_node_ids") && !isParamValid("secondary"))
    paramError("secondary", "Either secondary or secondary_node_ids must be provided.");
 
  // user provided conflicting params
  if (isParamValid("secondary_node_ids") && isParamValid("secondary"))
    paramError("secondary",
               "Both 'secondary' and 'secondary_node_ids' parameters are set. They are mutually "
               "exclusive.");
 
  // Fill in _connected_nodes, which defines local secondary nodes in the base class
  //
  // Note that later on when we pick the primary node from the secondary node set, we
  // will make sure to remove it from _connected_nodes
  _connected_nodes.clear();
 
  // user provided boundary name
  if (isParamValid("secondary"))
  {
    const auto & secondary_bnd = getParam<BoundaryName>("secondary");
    const auto & secondary_nodes = _mesh.getNodeList(_mesh.getBoundaryID(secondary_bnd));
 
    for (const auto & nid : secondary_nodes)
      if (_mesh.nodeRef(nid).processor_id() == _subproblem.processor_id())
        _connected_nodes.push_back(nid);
  }
  // user provided node ids
  else if (isParamValid("secondary_node_ids"))
  {
    const auto & secondary_node_ids = getParam<std::vector<unsigned int>>("secondary_node_ids");
    for (const auto & nid : secondary_node_ids)
      if (_mesh.queryNodePtr(nid) &&
          _mesh.nodeRef(nid).processor_id() == _subproblem.processor_id())
        _connected_nodes.push_back(nid);
  }
}
 
void
EqualValueBoundaryConstraint::pickPrimaryNode()
{
  // user provided nothing
  if (isParamValid("primary") && isParamValid("primary_node_coord"))
    mooseError(
        "Both 'primary' and 'primary_node_coord' parameters are set. They are mutually exclusive.");
 
  dof_id_type primary_node_id = Node::invalid_id;
 
  // user provided primary node coordinates
  if (isParamValid("primary_node_coord"))
    primary_node_id = getPrimaryNodeIDByCoord();
  // user provided primary node id
  else if (isParamValid("primary"))
    primary_node_id = getParam<unsigned int>("primary");
  // no primary node provided, so we pick one from secondary nodes
  else
  {
    // If the user provided secondary node ids, set primary node to first one
    if (isParamValid("secondary_node_ids"))
    {
      if (_connected_nodes.size())
        primary_node_id = _connected_nodes[0];
    }
    // otherwise, we pick the minimum node id from the secondary node set
    else
    {
      if (_connected_nodes.size())
        primary_node_id = (*std::min_element(_connected_nodes.begin(), _connected_nodes.end()));
      _mesh.comm().min(primary_node_id);
    }
  }
 
  mooseAssert(primary_node_id != Node::invalid_id, "We should have found a primary node");
 
  // Remove primary node from secondary nodes
  auto it = std::find(_connected_nodes.begin(), _connected_nodes.end(), primary_node_id);
  if (it != _connected_nodes.end())
    _connected_nodes.erase(it);
 
  // Store primary node id
  _primary_node_vector.resize(1);
  _primary_node_vector[0] = primary_node_id;
}
 
dof_id_type
EqualValueBoundaryConstraint::getPrimaryNodeIDByCoord() const
{
  const Real eps = libMesh::TOLERANCE;
  std::unordered_set<dof_id_type> local_primary_node_ids;
  const auto & primary_node_coord = getParam<Point>("primary_node_coord");
 
  // Gather local candidates
  for (const auto & bnd_node : *_mesh.getBoundaryNodeRange())
  {
    if ((*(bnd_node->_node) - primary_node_coord).norm() < eps)
    {
      // The primary node we found should belong to the secondary node set
      //
      // Note that we do not immediately break once a match is found because in theory, although it
      // is extremely unlikely, there could be multiple coinciding nodes on the secondary boundary
      // that match the provided coordinates. In that case, we would like to emit a meaning error
      // message.
      if (std::find(_connected_nodes.begin(), _connected_nodes.end(), bnd_node->_node->id()) !=
          _connected_nodes.end())
        local_primary_node_ids.insert(bnd_node->_node->id());
    }
  }
 
  // Gather all candidates from all ranks
  const std::vector<dof_id_type> local_node_vec(local_primary_node_ids.begin(),
                                                local_primary_node_ids.end());
  std::vector<std::vector<dof_id_type>> gathered_node_vecs;
  _mesh.comm().allgather(local_node_vec, gathered_node_vecs);
 
  // Deduplicate
  std::unordered_set<dof_id_type> global_primary_node_ids;
  for (const auto & vec : gathered_node_vecs)
    global_primary_node_ids.insert(vec.begin(), vec.end());
 
  // Make sure we found one and exactly one
  if (global_primary_node_ids.size() == 0)
    mooseError("Couldn't find a node ID for the specified primary_node_coord.");
  else if (global_primary_node_ids.size() > 1)
    mooseError("Multiple nodes found for the specified primary_node_coord.");
 
  return *global_primary_node_ids.begin();
}
 
void
EqualValueBoundaryConstraint::ghostPrimary()
{
  mooseAssert(_primary_node_vector.size() == 1, type() + " should have exactly one primary node");
 
  const auto primary_node_id = _primary_node_vector.front();
  const auto connected_elem_ids =
      gatherAndRetainConnectedElems(_fe_problem.mesh(), {primary_node_id});
  const auto primary_elem_id = connected_elem_ids.front();
  _subproblem.addGhostedElem(primary_elem_id);
 
  if (const auto displaced_problem = _fe_problem.getDisplacedProblem())
  {
    const auto displaced_connected_elem_ids =
        gatherAndRetainConnectedElems(displaced_problem->mesh(), {primary_node_id});
    const auto displaced_primary_elem_id = displaced_connected_elem_ids.front();
    if (displaced_primary_elem_id != primary_elem_id)
      mooseError("Reference and displaced meshes selected different primary elements");
  }
}
 
Real
EqualValueBoundaryConstraint::computeQpResidual(Moose::ConstraintType type)
{
  switch (type)
  {
    case Moose::Secondary:
      return (_u_secondary[_i] - _u_primary[_j]) * _penalty;
    case Moose::Primary:
      return (_u_primary[_j] - _u_secondary[_i]) * _penalty;
  }
  return 0.;
}
 
Real
EqualValueBoundaryConstraint::computeQpJacobian(Moose::ConstraintJacobianType type)
{
  switch (type)
  {
    case Moose::SecondarySecondary:
      return _penalty;
    case Moose::SecondaryPrimary:
      return -_penalty;
    case Moose::PrimaryPrimary:
      return _penalty;
    case Moose::PrimarySecondary:
      return -_penalty;
    default:
      mooseError("Unsupported type");
      break;
  }
  return 0.;
}