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             : // C++ includes
      20             : #include <limits.h>
      21             : 
      22             : namespace // Anonymous namespace for helpers
      23             : {
      24             : /**
      25             :  * Specific weak ordering for Elem *'s to be used in a set.
      26             :  * We use the id, but first sort by level.  This guarantees
      27             :  * when traversing the set from beginning to end the lower
      28             :  * level (parent) elements are encountered first.
      29             :  *
      30             :  * This was swiped from libMesh mesh_communication.C, and ought to be
      31             :  * replaced with libMesh::CompareElemIdsByLevel just as soon as I refactor to
      32             :  * create that - @roystgnr
      33             :  */
      34             : struct CompareElemsByLevel
      35             : {
      36           6 :   bool operator()(const Elem * a, const Elem * b) const
      37             :   {
      38             :     libmesh_assert(a);
      39             :     libmesh_assert(b);
      40           6 :     const unsigned int al = a->level(), bl = b->level();
      41           6 :     const dof_id_type aid = a->id(), bid = b->id();
      42             : 
      43           6 :     return (al == bl) ? aid < bid : al < bl;
      44             :   }
      45             : };
      46             : 
      47             : } // anonymous namespace
      48             : 
      49             : registerMooseObject("MooseApp", EqualValueBoundaryConstraint);
      50             : 
      51             : InputParameters
      52       14309 : EqualValueBoundaryConstraint::validParams()
      53             : {
      54       14309 :   InputParameters params = NodalConstraint::validParams();
      55       14309 :   params.addClassDescription(
      56             :       "Constraint for enforcing that variables on each side of a boundary are equivalent.");
      57       42927 :   params.addParam<unsigned int>(
      58             :       "primary",
      59       28618 :       std::numeric_limits<unsigned int>::max(),
      60             :       "The ID of the primary node. If no ID is provided, first node of secondary set is chosen.");
      61       14309 :   params.addParam<std::vector<unsigned int>>(
      62             :       "secondary_node_ids", {}, "The IDs of the secondary node");
      63       14309 :   params.addParam<BoundaryName>(
      64             :       "secondary", "NaN", "The boundary ID associated with the secondary side");
      65       14309 :   params.addRequiredParam<Real>("penalty", "The penalty used for the boundary term");
      66       14309 :   return params;
      67           0 : }
      68             : 
      69          22 : EqualValueBoundaryConstraint::EqualValueBoundaryConstraint(const InputParameters & parameters)
      70             :   : NodalConstraint(parameters),
      71          22 :     _primary_node_id(getParam<unsigned int>("primary")),
      72          22 :     _secondary_node_ids(getParam<std::vector<unsigned int>>("secondary_node_ids")),
      73          22 :     _secondary_node_set_id(getParam<BoundaryName>("secondary")),
      74          44 :     _penalty(getParam<Real>("penalty"))
      75             : {
      76          22 :   updateConstrainedNodes();
      77          22 : }
      78             : 
      79             : void
      80          18 : EqualValueBoundaryConstraint::meshChanged()
      81             : {
      82          18 :   updateConstrainedNodes();
      83          18 : }
      84             : 
      85             : void
      86          40 : EqualValueBoundaryConstraint::updateConstrainedNodes()
      87             : {
      88          40 :   _primary_node_vector.clear();
      89          40 :   _connected_nodes.clear();
      90             : 
      91          40 :   if ((_secondary_node_ids.size() == 0) && (_secondary_node_set_id == "NaN"))
      92           0 :     mooseError("Please specify secondary node ids or boundary id.");
      93          40 :   else if ((_secondary_node_ids.size() == 0) && (_secondary_node_set_id != "NaN"))
      94             :   {
      95             :     std::vector<dof_id_type> nodelist =
      96          40 :         _mesh.getNodeList(_mesh.getBoundaryID(_secondary_node_set_id));
      97          40 :     std::vector<dof_id_type>::iterator in;
      98             : 
      99             :     // Set primary node to first node of the secondary node set if no primary node id is provided
     100             :     //_primary_node_vector defines primary nodes in the base class
     101          40 :     if (_primary_node_id == std::numeric_limits<unsigned int>::max())
     102             :     {
     103           0 :       in = std::min_element(nodelist.begin(), nodelist.end());
     104           0 :       dof_id_type node_id = (in == nodelist.end()) ? DofObject::invalid_id : *in;
     105           0 :       _communicator.min(node_id);
     106           0 :       _primary_node_vector.push_back(node_id);
     107             :     }
     108             :     else
     109          40 :       _primary_node_vector.push_back(_primary_node_id);
     110             : 
     111             :     // Fill in _connected_nodes, which defines secondary nodes in the base class
     112         532 :     for (in = nodelist.begin(); in != nodelist.end(); ++in)
     113             :     {
     114         944 :       if ((*in != _primary_node_vector[0]) &&
     115         452 :           (_mesh.nodeRef(*in).processor_id() == _subproblem.processor_id()))
     116         354 :         _connected_nodes.push_back(*in);
     117             :     }
     118          40 :   }
     119           0 :   else if ((_secondary_node_ids.size() != 0) && (_secondary_node_set_id == "NaN"))
     120             :   {
     121           0 :     if (_primary_node_id == std::numeric_limits<unsigned int>::max())
     122           0 :       _primary_node_vector.push_back(
     123           0 :           _secondary_node_ids[0]); //_primary_node_vector defines primary nodes in the base class
     124             : 
     125             :     // Fill in _connected_nodes, which defines secondary nodes in the base class
     126           0 :     for (const auto & dof : _secondary_node_ids)
     127             :     {
     128           0 :       if (_mesh.queryNodePtr(dof) &&
     129           0 :           (_mesh.nodeRef(dof).processor_id() == _subproblem.processor_id()) &&
     130           0 :           (dof != _primary_node_vector[0]))
     131           0 :         _connected_nodes.push_back(dof);
     132             :     }
     133             :   }
     134             : 
     135          40 :   const auto & node_to_elem_map = _mesh.nodeToElemMap();
     136          40 :   auto node_to_elem_pair = node_to_elem_map.find(_primary_node_vector[0]);
     137             : 
     138          40 :   bool found_elems = (node_to_elem_pair != node_to_elem_map.end());
     139             : 
     140             :   // Add elements connected to primary node to Ghosted Elements.
     141             : 
     142             :   // On a distributed mesh, these elements might have already been
     143             :   // remoted, in which case we need to gather them back first.
     144          40 :   if (!_mesh.getMesh().is_serial())
     145             :   {
     146             : #ifndef NDEBUG
     147             :     bool someone_found_elems = found_elems;
     148             :     _mesh.getMesh().comm().max(someone_found_elems);
     149             :     mooseAssert(someone_found_elems, "Missing entry in node to elem map");
     150             : #endif
     151             : 
     152           2 :     std::set<Elem *, CompareElemsByLevel> primary_elems_to_ghost;
     153           2 :     std::set<Node *> nodes_to_ghost;
     154           2 :     if (found_elems)
     155             :     {
     156           6 :       for (dof_id_type id : node_to_elem_pair->second)
     157             :       {
     158           4 :         Elem * elem = _mesh.queryElemPtr(id);
     159           4 :         if (elem)
     160             :         {
     161           4 :           primary_elems_to_ghost.insert(elem);
     162             : 
     163           4 :           const unsigned int n_nodes = elem->n_nodes();
     164          20 :           for (unsigned int n = 0; n != n_nodes; ++n)
     165          16 :             nodes_to_ghost.insert(elem->node_ptr(n));
     166             :         }
     167             :       }
     168             :     }
     169             : 
     170             :     // Send nodes first since elements need them
     171           2 :     _mesh.getMesh().comm().allgather_packed_range(&_mesh.getMesh(),
     172             :                                                   nodes_to_ghost.begin(),
     173             :                                                   nodes_to_ghost.end(),
     174             :                                                   libMesh::null_output_iterator<Node>());
     175             : 
     176           2 :     _mesh.getMesh().comm().allgather_packed_range(&_mesh.getMesh(),
     177             :                                                   primary_elems_to_ghost.begin(),
     178             :                                                   primary_elems_to_ghost.end(),
     179             :                                                   libMesh::null_output_iterator<Elem>());
     180             : 
     181           2 :     _mesh.update(); // Rebuild node_to_elem_map
     182             : 
     183             :     // Find elems again now that we know they're there
     184           2 :     const auto & new_node_to_elem_map = _mesh.nodeToElemMap();
     185           2 :     node_to_elem_pair = new_node_to_elem_map.find(_primary_node_vector[0]);
     186           2 :     found_elems = (node_to_elem_pair != new_node_to_elem_map.end());
     187           2 :   }
     188             : 
     189          40 :   if (!found_elems)
     190           0 :     mooseError("Couldn't find any elements connected to primary node");
     191             : 
     192          40 :   const std::vector<dof_id_type> & elems = node_to_elem_pair->second;
     193             : 
     194          40 :   if (elems.size() == 0)
     195           0 :     mooseError("Couldn't find any elements connected to primary node");
     196          40 :   _subproblem.addGhostedElem(elems[0]);
     197          40 : }
     198             : 
     199             : Real
     200        3204 : EqualValueBoundaryConstraint::computeQpResidual(Moose::ConstraintType type)
     201             : {
     202        3204 :   switch (type)
     203             :   {
     204        1602 :     case Moose::Secondary:
     205        1602 :       return (_u_secondary[_i] - _u_primary[_j]) * _penalty;
     206        1602 :     case Moose::Primary:
     207        1602 :       return (_u_primary[_j] - _u_secondary[_i]) * _penalty;
     208             :   }
     209           0 :   return 0.;
     210             : }
     211             : 
     212             : Real
     213        1752 : EqualValueBoundaryConstraint::computeQpJacobian(Moose::ConstraintJacobianType type)
     214             : {
     215        1752 :   switch (type)
     216             :   {
     217         438 :     case Moose::SecondarySecondary:
     218         438 :       return _penalty;
     219         438 :     case Moose::SecondaryPrimary:
     220         438 :       return -_penalty;
     221         438 :     case Moose::PrimaryPrimary:
     222         438 :       return _penalty;
     223         438 :     case Moose::PrimarySecondary:
     224         438 :       return -_penalty;
     225           0 :     default:
     226           0 :       mooseError("Unsupported type");
     227             :       break;
     228             :   }
     229             :   return 0.;
     230             : }