Line data    Source code

       1             : //* This file is part of the MOOSE framework
       2             : //* https://mooseframework.inl.gov
       3             : //*
       4             : //* All rights reserved, see COPYRIGHT for full restrictions
       5             : //* https://github.com/idaholab/moose/blob/master/COPYRIGHT
       6             : //*
       7             : //* Licensed under LGPL 2.1, please see LICENSE for details
       8             : //* https://www.gnu.org/licenses/lgpl-2.1.html
       9             : 
      10             : // MOOSE includes
      11             : #include "EqualValueBoundaryConstraint.h"
      12             : #include "MooseMesh.h"
      13             : 
      14             : #include "libmesh/null_output_iterator.h"
      15             : #include "libmesh/parallel.h"
      16             : #include "libmesh/parallel_elem.h"
      17             : #include "libmesh/parallel_node.h"
      18             : 
      19             : namespace // Anonymous namespace for helpers
      20             : {
      21             : /**
      22             :  * Specific weak ordering for Elem *'s to be used in a set.
      23             :  * We use the id, but first sort by level.  This guarantees
      24             :  * when traversing the set from beginning to end the lower
      25             :  * level (parent) elements are encountered first.
      26             :  *
      27             :  * This was swiped from libMesh mesh_communication.C, and ought to be
      28             :  * replaced with libMesh::CompareElemIdsByLevel just as soon as I refactor to
      29             :  * create that - @roystgnr
      30             :  */
      31             : struct CompareElemsByLevel
      32             : {
      33           9 :   bool operator()(const Elem * a, const Elem * b) const
      34             :   {
      35             :     libmesh_assert(a);
      36             :     libmesh_assert(b);
      37           9 :     const unsigned int al = a->level(), bl = b->level();
      38           9 :     const dof_id_type aid = a->id(), bid = b->id();
      39             : 
      40           9 :     return (al == bl) ? aid < bid : al < bl;
      41             :   }
      42             : };
      43             : 
      44             : } // anonymous namespace
      45             : 
      46             : registerMooseObject("MooseApp", EqualValueBoundaryConstraint);
      47             : 
      48             : InputParameters
      49        3138 : EqualValueBoundaryConstraint::validParams()
      50             : {
      51        3138 :   InputParameters params = NodalConstraint::validParams();
      52        6276 :   params.addClassDescription(
      53             :       "Constraint for enforcing that variables on each side of a boundary are equivalent.");
      54       12552 :   params.addParam<unsigned int>(
      55             :       "primary",
      56             :       "The ID of the primary node. If no ID is provided, first node of secondary set is chosen.");
      57       12552 :   params.addParam<Point>("primary_node_coord", "Coordinates of the primary node to locate.");
      58       12552 :   params.addParam<std::vector<unsigned int>>("secondary_node_ids", "The IDs of the secondary node");
      59       12552 :   params.addParam<BoundaryName>("secondary", "The boundary ID associated with the secondary side");
      60        9414 :   params.addRequiredParam<Real>("penalty", "The penalty used for the boundary term");
      61        3138 :   return params;
      62           0 : }
      63             : 
      64          43 : EqualValueBoundaryConstraint::EqualValueBoundaryConstraint(const InputParameters & parameters)
      65          86 :   : NodalConstraint(parameters), _penalty(getParam<Real>("penalty"))
      66             : {
      67          43 :   updateConstrainedNodes();
      68          34 : }
      69             : 
      70             : void
      71          18 : EqualValueBoundaryConstraint::meshChanged()
      72             : {
      73          18 :   updateConstrainedNodes();
      74          18 : }
      75             : 
      76             : void
      77          61 : EqualValueBoundaryConstraint::updateConstrainedNodes()
      78             : {
      79          61 :   populateSecondaryNodes();
      80          61 :   pickPrimaryNode();
      81          52 :   ghostPrimary();
      82          52 : }
      83             : 
      84             : void
      85          61 : EqualValueBoundaryConstraint::populateSecondaryNodes()
      86             : {
      87             :   // user provided nothing
      88         360 :   if (!isParamValid("secondary_node_ids") && !isParamValid("secondary"))
      89           0 :     paramError("secondary", "Either secondary or secondary_node_ids must be provided.");
      90             : 
      91             :   // user provided conflicting params
      92         189 :   if (isParamValid("secondary_node_ids") && isParamValid("secondary"))
      93           0 :     paramError("secondary",
      94             :                "Both 'secondary' and 'secondary_node_ids' parameters are set. They are mutually "
      95             :                "exclusive.");
      96             : 
      97             :   // Fill in _connected_nodes, which defines local secondary nodes in the base class
      98             :   //
      99             :   // Note that later on when we pick the primary node from the secondary node set, we
     100             :   // will make sure to remove it from _connected_nodes
     101          61 :   _connected_nodes.clear();
     102             : 
     103             :   // user provided boundary name
     104         183 :   if (isParamValid("secondary"))
     105             :   {
     106         116 :     const auto & secondary_bnd = getParam<BoundaryName>("secondary");
     107          58 :     const auto & secondary_nodes = _mesh.getNodeList(_mesh.getBoundaryID(secondary_bnd));
     108             : 
     109         672 :     for (const auto & nid : secondary_nodes)
     110         614 :       if (_mesh.nodeRef(nid).processor_id() == _subproblem.processor_id())
     111         490 :         _connected_nodes.push_back(nid);
     112             :   }
     113             :   // user provided node ids
     114           9 :   else if (isParamValid("secondary_node_ids"))
     115             :   {
     116           6 :     const auto & secondary_node_ids = getParam<std::vector<unsigned int>>("secondary_node_ids");
     117          21 :     for (const auto & nid : secondary_node_ids)
     118          33 :       if (_mesh.queryNodePtr(nid) &&
     119          15 :           _mesh.nodeRef(nid).processor_id() == _subproblem.processor_id())
     120          15 :         _connected_nodes.push_back(nid);
     121             :   }
     122          61 : }
     123             : 
     124             : void
     125          61 : EqualValueBoundaryConstraint::pickPrimaryNode()
     126             : {
     127             :   // user provided nothing
     128         269 :   if (isParamValid("primary") && isParamValid("primary_node_coord"))
     129           3 :     mooseError(
     130             :         "Both 'primary' and 'primary_node_coord' parameters are set. They are mutually exclusive.");
     131             : 
     132          58 :   dof_id_type primary_node_id = Node::invalid_id;
     133             : 
     134             :   // user provided primary node coordinates
     135         174 :   if (isParamValid("primary_node_coord"))
     136          18 :     primary_node_id = getPrimaryNodeIDByCoord();
     137             :   // user provided primary node id
     138         120 :   else if (isParamValid("primary"))
     139         120 :     primary_node_id = getParam<unsigned int>("primary");
     140             :   // no primary node provided, so we pick one from secondary nodes
     141             :   else
     142             :   {
     143             :     // If the user provided secondary node ids, set primary node to first one
     144           0 :     if (isParamValid("secondary_node_ids"))
     145             :     {
     146           0 :       if (_connected_nodes.size())
     147           0 :         primary_node_id = _connected_nodes[0];
     148             :     }
     149             :     // otherwise, we pick the minimum node id from the secondary node set
     150             :     else
     151             :     {
     152           0 :       if (_connected_nodes.size())
     153           0 :         primary_node_id = (*std::min_element(_connected_nodes.begin(), _connected_nodes.end()));
     154           0 :       _mesh.comm().min(primary_node_id);
     155             :     }
     156             :   }
     157             : 
     158             :   mooseAssert(primary_node_id != Node::invalid_id, "We should have found a primary node");
     159             : 
     160             :   // Remove primary node from secondary nodes
     161          52 :   auto it = std::find(_connected_nodes.begin(), _connected_nodes.end(), primary_node_id);
     162          52 :   if (it != _connected_nodes.end())
     163          40 :     _connected_nodes.erase(it);
     164             : 
     165             :   // Store primary node id
     166          52 :   _primary_node_vector.resize(1);
     167          52 :   _primary_node_vector[0] = primary_node_id;
     168          52 : }
     169             : 
     170             : dof_id_type
     171          18 : EqualValueBoundaryConstraint::getPrimaryNodeIDByCoord() const
     172             : {
     173          18 :   const Real eps = libMesh::TOLERANCE;
     174          18 :   std::unordered_set<dof_id_type> local_primary_node_ids;
     175          36 :   const auto & primary_node_coord = getParam<Point>("primary_node_coord");
     176             : 
     177             :   // Gather local candidates
     178         476 :   for (const auto & bnd_node : *_mesh.getBoundaryNodeRange())
     179             :   {
     180         458 :     if ((*(bnd_node->_node) - primary_node_coord).norm() < eps)
     181             :     {
     182             :       // The primary node we found should belong to the secondary node set
     183             :       //
     184             :       // Note that we do not immediately break once a match is found because in theory, although it
     185             :       // is extremely unlikely, there could be multiple coinciding nodes on the secondary boundary
     186             :       // that match the provided coordinates. In that case, we would like to emit a meaning error
     187             :       // message.
     188          30 :       if (std::find(_connected_nodes.begin(), _connected_nodes.end(), bnd_node->_node->id()) !=
     189          60 :           _connected_nodes.end())
     190          27 :         local_primary_node_ids.insert(bnd_node->_node->id());
     191             :     }
     192             :   }
     193             : 
     194             :   // Gather all candidates from all ranks
     195             :   const std::vector<dof_id_type> local_node_vec(local_primary_node_ids.begin(),
     196          18 :                                                 local_primary_node_ids.end());
     197          18 :   std::vector<std::vector<dof_id_type>> gathered_node_vecs;
     198          18 :   _mesh.comm().allgather(local_node_vec, gathered_node_vecs);
     199             : 
     200             :   // Deduplicate
     201          18 :   std::unordered_set<dof_id_type> global_primary_node_ids;
     202          42 :   for (const auto & vec : gathered_node_vecs)
     203          24 :     global_primary_node_ids.insert(vec.begin(), vec.end());
     204             : 
     205             :   // Make sure we found one and exactly one
     206          18 :   if (global_primary_node_ids.size() == 0)
     207           3 :     mooseError("Couldn't find a node ID for the specified primary_node_coord.");
     208          15 :   else if (global_primary_node_ids.size() > 1)
     209           3 :     mooseError("Multiple nodes found for the specified primary_node_coord.");
     210             : 
     211          24 :   return *global_primary_node_ids.begin();
     212          12 : }
     213             : 
     214             : void
     215          52 : EqualValueBoundaryConstraint::ghostPrimary()
     216             : {
     217          52 :   const auto & node_to_elem_map = _mesh.nodeToElemMap();
     218          52 :   auto node_to_elem_pair = node_to_elem_map.find(_primary_node_vector[0]);
     219          52 :   bool found_elems = (node_to_elem_pair != node_to_elem_map.end());
     220             : 
     221             :   // Add elements connected to primary node to Ghosted Elements.
     222             : 
     223             :   // On a distributed mesh, these elements might have already been
     224             :   // remoted, in which case we need to gather them back first.
     225          52 :   if (!_mesh.getMesh().is_serial())
     226             :   {
     227             : #ifndef NDEBUG
     228             :     bool someone_found_elems = found_elems;
     229             :     _mesh.getMesh().comm().max(someone_found_elems);
     230             :     mooseAssert(someone_found_elems, "Missing entry in node to elem map");
     231             : #endif
     232             : 
     233           4 :     std::set<Elem *, CompareElemsByLevel> primary_elems_to_ghost;
     234           4 :     std::set<Node *> nodes_to_ghost;
     235           4 :     if (found_elems)
     236             :     {
     237          11 :       for (dof_id_type id : node_to_elem_pair->second)
     238             :       {
     239           7 :         Elem * elem = _mesh.queryElemPtr(id);
     240           7 :         if (elem)
     241             :         {
     242           7 :           primary_elems_to_ghost.insert(elem);
     243             : 
     244           7 :           const unsigned int n_nodes = elem->n_nodes();
     245          35 :           for (unsigned int n = 0; n != n_nodes; ++n)
     246          28 :             nodes_to_ghost.insert(elem->node_ptr(n));
     247             :         }
     248             :       }
     249             :     }
     250             : 
     251             :     // Send nodes first since elements need them
     252           4 :     _mesh.getMesh().comm().allgather_packed_range(&_mesh.getMesh(),
     253             :                                                   nodes_to_ghost.begin(),
     254             :                                                   nodes_to_ghost.end(),
     255             :                                                   libMesh::null_output_iterator<Node>());
     256             : 
     257           4 :     _mesh.getMesh().comm().allgather_packed_range(&_mesh.getMesh(),
     258             :                                                   primary_elems_to_ghost.begin(),
     259             :                                                   primary_elems_to_ghost.end(),
     260             :                                                   libMesh::null_output_iterator<Elem>());
     261             : 
     262             :     // After allgather_packed_range(), rebuild internal connectivity.
     263             :     // This updates ghost nodes/elements across processors and reconstructs node_to_elem_map.
     264           4 :     _mesh.update();
     265             : 
     266             :     // Find elems again now that we know they're there
     267           4 :     const auto & new_node_to_elem_map = _mesh.nodeToElemMap();
     268           4 :     node_to_elem_pair = new_node_to_elem_map.find(_primary_node_vector[0]);
     269           4 :     found_elems = (node_to_elem_pair != new_node_to_elem_map.end());
     270           4 :   }
     271             : 
     272          52 :   if (!found_elems)
     273           0 :     mooseError("Couldn't find any elements connected to primary node");
     274             : 
     275          52 :   const std::vector<dof_id_type> & elems = node_to_elem_pair->second;
     276             : 
     277          52 :   if (elems.size() == 0)
     278           0 :     mooseError("Couldn't find any elements connected to primary node");
     279          52 :   _subproblem.addGhostedElem(elems[0]);
     280          52 : }
     281             : 
     282             : Real
     283        5532 : EqualValueBoundaryConstraint::computeQpResidual(Moose::ConstraintType type)
     284             : {
     285        5532 :   switch (type)
     286             :   {
     287        2766 :     case Moose::Secondary:
     288        2766 :       return (_u_secondary[_i] - _u_primary[_j]) * _penalty;
     289        2766 :     case Moose::Primary:
     290        2766 :       return (_u_primary[_j] - _u_secondary[_i]) * _penalty;
     291             :   }
     292           0 :   return 0.;
     293             : }
     294             : 
     295             : Real
     296        2376 : EqualValueBoundaryConstraint::computeQpJacobian(Moose::ConstraintJacobianType type)
     297             : {
     298        2376 :   switch (type)
     299             :   {
     300         594 :     case Moose::SecondarySecondary:
     301         594 :       return _penalty;
     302         594 :     case Moose::SecondaryPrimary:
     303         594 :       return -_penalty;
     304         594 :     case Moose::PrimaryPrimary:
     305         594 :       return _penalty;
     306         594 :     case Moose::PrimarySecondary:
     307         594 :       return -_penalty;
     308           0 :     default:
     309           0 :       mooseError("Unsupported type");
     310             :       break;
     311             :   }
     312             :   return 0.;
     313             : }