Demonstrate bound constraint enforcement using nodal kernels

Refs #2999

framework/include/nodalkernels/CoupledForceNodalKernel.h

@@ -0,0 +1,42 @@
+//* This file is part of the MOOSE framework
+//* https://www.mooseframework.org
+//*
+//* 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
+
+#pragma once
+
+#include "NodalKernel.h"
+
+// Forward Declarations
+class CoupledForceNodalKernel;
+
+template <>
+InputParameters validParams<CoupledForceNodalKernel>();
+
+/**
+ * Adds a force proportional to the value of the coupled variable
+ */
+class CoupledForceNodalKernel : public NodalKernel
+{
+public:
+  CoupledForceNodalKernel(const InputParameters & parameters);
+
+protected:
+  virtual Real computeQpResidual() override;
+  virtual Real computeQpJacobian() override;
+  virtual Real computeQpOffDiagJacobian(unsigned int jvar) override;
+
+private:
+  /// The number of the coupled variable
+  const unsigned int _v_var;
+
+  /// The value of the coupled variable
+  const VariableValue & _v;
+
+  /// A multiplicative factor for computing the coupled force
+  const Real _coef;
+};

framework/include/nodalkernels/LowerBoundNodalKernel.h

@@ -0,0 +1,45 @@
+//* This file is part of the MOOSE framework
+//* https://www.mooseframework.org
+//*
+//* 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
+
+#pragma once
+
+#include "NodalKernel.h"
+
+// Forward Declarations
+class LowerBoundNodalKernel;
+
+template <>
+InputParameters validParams<LowerBoundNodalKernel>();
+
+/**
+ * Class used to enforce a lower bound on a coupled variable
+ */
+class LowerBoundNodalKernel : public NodalKernel
+{
+public:
+  LowerBoundNodalKernel(const InputParameters & parameters);
+
+protected:
+  virtual Real computeQpResidual() override;
+  virtual Real computeQpJacobian() override;
+  virtual Real computeQpOffDiagJacobian(unsigned int jvar) override;
+
+private:
+  /// The number of the coupled variable
+  const unsigned int _v_var;
+
+  /// The value of the coupled variable
+  const VariableValue & _v;
+
+  /// Boundaries on which we should not execute this object
+  std::set<BoundaryID> _bnd_ids;
+
+  /// The lower bound on the coupled variable
+  const Real _lower_bound;
+};

framework/include/postprocessors/LMActiveSetSize.h

@@ -0,0 +1,56 @@
+//* This file is part of the MOOSE framework
+//* https://www.mooseframework.org
+//*
+//* 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
+
+#pragma once
+
+#include "GeneralPostprocessor.h"
+
+// Forward Declarations
+class LMActiveSetSize;
+class MooseVariableFEBase;
+namespace libMesh
+{
+class MeshBase;
+}
+
+template <>
+InputParameters validParams<LMActiveSetSize>();
+
+class LMActiveSetSize : public GeneralPostprocessor
+{
+public:
+  LMActiveSetSize(const InputParameters & parameters);
+
+  void initialize() override;
+  void execute() override;
+
+  PostprocessorValue getValue() override;
+
+private:
+  /// MOOSE variable we compute the contact set from
+  const MooseVariableFEBase & _var;
+
+  /// The libmesh mesh
+  const MeshBase & _mesh;
+
+  /// Whether we are subdomain restricting the active set search
+  const bool _subdomain_restricted;
+
+  /// An optional subdomain over which to query degrees of freedom
+  const SubdomainID _subdomain_id;
+
+  /// The tolerance used to decide whether the variable indicates contact
+  const Real _value;
+
+  /// Represents the number of values in contact
+  unsigned int _count;
+
+  /// The comparison to perform
+  const MooseEnum _comparator;
+};

framework/src/nodalkernels/CoupledForceNodalKernel.C

@@ -0,0 +1,57 @@
+//* This file is part of the MOOSE framework
+//* https://www.mooseframework.org
+//*
+//* 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
+
+#include "CoupledForceNodalKernel.h"
+
+registerMooseObject("MooseApp", CoupledForceNodalKernel);
+
+template <>
+InputParameters
+validParams<CoupledForceNodalKernel>()
+{
+  InputParameters params = validParams<NodalKernel>();
+  params.addClassDescription("Adds a force proportional to the value of the coupled variable");
+  params.addRequiredCoupledVar("v", "The coupled variable which provides the force");
+  params.addParam<Real>(
+      "coef", 1.0, "Coefficent ($\\sigma$) multiplier for the coupled force term.");
+
+  return params;
+}
+
+CoupledForceNodalKernel::CoupledForceNodalKernel(const InputParameters & parameters)
+  : NodalKernel(parameters),
+    _v_var(coupled("v")),
+    _v(coupledValue("v")),
+    _coef(getParam<Real>("coef"))
+{
+  if (_var.number() == _v_var)
+    mooseError(
+        "Coupled variable 'v' needs to be different from 'variable' with CoupledForceNodalKernel, "
+        "consider using Reaction or somethig similar");
+}
+
+Real
+CoupledForceNodalKernel::computeQpResidual()
+{
+  return -_coef * _v[_qp];
+}
+
+Real
+CoupledForceNodalKernel::computeQpJacobian()
+{
+  return 0;
+}
+
+Real
+CoupledForceNodalKernel::computeQpOffDiagJacobian(unsigned int jvar)
+{
+  if (jvar == _v_var)
+    return -_coef;
+  return 0.0;
+}

framework/src/nodalkernels/LowerBoundNodalKernel.C

@@ -0,0 +1,80 @@
+//* This file is part of the MOOSE framework
+//* https://www.mooseframework.org
+//*
+//* 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
+
+#include "LowerBoundNodalKernel.h"
+
+registerMooseObject("MooseApp", LowerBoundNodalKernel);
+
+template <>
+InputParameters
+validParams<LowerBoundNodalKernel>()
+{
+  InputParameters params = validParams<NodalKernel>();
+  params.addClassDescription("Used to prevent a coupled variable from going below a lower bound");
+  params.addRequiredCoupledVar(
+      "v", "The coupled variable we require to be greater than the lower bound");
+  params.addParam<Real>("lower_bound", 0, "The lower bound on the coupled variable");
+  params.addParam<std::vector<BoundaryName>>(
+      "exclude_boundaries",
+      "Boundaries on which not to execute the NodalKernel. This can be useful for avoiding "
+      "singuarility in the matrix in case a constraint is active in the same place that a "
+      "DirichletBC is set");
+  return params;
+}
+
+LowerBoundNodalKernel::LowerBoundNodalKernel(const InputParameters & parameters)
+  : NodalKernel(parameters),
+    _v_var(coupled("v")),
+    _v(coupledValue("v")),
+    _lower_bound(getParam<Real>("lower_bound"))
+{
+  if (_var.number() == _v_var)
+    mooseError("Coupled variable 'v' needs to be different from 'variable' with "
+               "LowerBoundNodalKernel");
+
+  const auto & bnd_names = getParam<std::vector<BoundaryName>>("exclude_boundaries");
+  for (const auto & bnd_name : bnd_names)
+    _bnd_ids.insert(_mesh.getBoundaryID(bnd_name));
+}
+
+Real
+LowerBoundNodalKernel::computeQpResidual()
+{
+  for (auto bnd_id : _bnd_ids)
+    if (_mesh.isBoundaryNode(_current_node->id(), bnd_id))
+      return _u[_qp];
+
+  return std::min(_u[_qp], _v[_qp] - _lower_bound);
+}
+
+Real
+LowerBoundNodalKernel::computeQpJacobian()
+{
+  for (auto bnd_id : _bnd_ids)
+    if (_mesh.isBoundaryNode(_current_node->id(), bnd_id))
+      return 1;
+
+  if (_u[_qp] <= _v[_qp] - _lower_bound)
+    return 1;
+  return 0;
+}
+
+Real
+LowerBoundNodalKernel::computeQpOffDiagJacobian(unsigned int jvar)
+{
+  for (auto bnd_id : _bnd_ids)
+    if (_mesh.isBoundaryNode(_current_node->id(), bnd_id))
+      return 0;
+
+  if (jvar == _v_var)
+    if (_v[_qp] - _lower_bound < _u[_qp])
+      return 1;
+
+  return 0.0;
+}

framework/src/nodalkernels/NodalKernel.C

@@ -187,19 +187,23 @@ NodalKernel::computeJacobian()
 void
 NodalKernel::computeOffDiagJacobian(unsigned int jvar)
 {
-  if (jvar == _var.number())
-    computeJacobian();
-  else
+  if (_var.isNodalDefined())
   {
-    _qp = 0;
-    Real cached_val = computeQpOffDiagJacobian(jvar);
-    dof_id_type cached_row = _var.nodalDofIndex();
-    // Note: this only works for Lagrange variables...
-    dof_id_type cached_col = _current_node->dof_number(_sys.number(), jvar, 0);
+    if (jvar == _var.number())
+      computeJacobian();
+    else
+    {
+      _qp = 0;
+      Real cached_val = computeQpOffDiagJacobian(jvar);
+      dof_id_type cached_row = _var.nodalDofIndex();
 
-    cached_val *= _var.scalingFactor();
+      // Note: this only works for equal order Lagrange variables...
+      dof_id_type cached_col = _current_node->dof_number(_sys.number(), jvar, 0);
+
+      cached_val *= _var.scalingFactor();
 
-    _assembly.cacheJacobianContribution(cached_row, cached_col, cached_val, _matrix_tags);
+      _assembly.cacheJacobianContribution(cached_row, cached_col, cached_val, _matrix_tags);
+    }
   }
 }