Support aux kernel evaluation of lower-d vars

framework/doc/content/syntax/AuxKernels/index.md

@@ -67,12 +67,12 @@ is done to allow the syntax to be consistent regardless of the AuxKernel flavor:
 
 In order to compute properties in the mortar sense, it is necessary to loop over the mortar segment
 mesh to spatially integrate variables. `MortarNodalAuxKernel`s offer this functionality where these "weighted" variables,
-which intervene in the computation of contact constraints and their residuals, can be coupled to generate the desired output value. 
-Therefore, if postprocessing of mortar quantities is required, nodal mortar auxiliary kernels can be employed. 
+which intervene in the computation of contact constraints and their residuals, can be coupled to generate the desired output value.
+Therefore, if postprocessing of mortar quantities is required, nodal mortar auxiliary kernels can be employed.
 Objects inheriting from `MortarNodalAuxKernel` allow for said operations on the mortar lower-dimensional domains featuring similar
 functionality to other nodal auxiliary kernels, including the possibility of computing quantities in an
 `incremental` manner.
- 
+
 ## Execute Flags
 
 AuxKernel objects inherit from the [SetupInterface.md] so they include the "execute_on" variable.
@@ -88,6 +88,16 @@ execution of a simulation. In this case, the `EXEC_LINEAR` flag should be remove
 `EXEC_INITIAL` flag should be added to perform the auxiliary variable calculation during the initial
 setup phase as well.
 
+## Populating lower-dimensional auxiliary variables
+
+Lower-dimensional auxiliary variables may be populated using boundary restricted
+auxiliary kernels. The boundary restriction of the aux kernel should be
+coincident with (a subset of) the lower-dimensional blocks that the
+lower-dimensional variable lives on. Using a boundary restricted auxiliary
+kernel as opposed to a lower-d block-restricted auxiliary kernel allows pulling
+in coincident face evaluations of higher-dimensional variables and material
+properties as well as evaluations of coupled lower-dimensional variables.
+
 ## Example A: Post processing with AuxKernel
 
 The following example is extracted from step 4 of the

framework/include/auxkernels/AuxKernel.h

@@ -116,6 +116,11 @@ class AuxKernelTempl : public MooseObject,
   template <typename T>
   const MaterialProperty<T> & getMaterialPropertyOlder(const std::string & name);
 
+  /**
+   * Insert the just computed values into the auxiliary solution vector
+   */
+  void insert();
+
 protected:
   /**
    * Compute and return the value of the aux variable.
@@ -214,6 +219,13 @@ class AuxKernelTempl : public MooseObject,
   /// reference to the solution vector of auxiliary system
   NumericVector<Number> & _solution;
 
+  /// The current lower dimensional element
+  const Elem * const & _current_lower_d_elem;
+
+  /// Whether we are computing for a lower dimensional variable using boundary restriction, e.g. a
+  /// variable whose block restriction is coincident with a higher-dimensional boundary face
+  const bool _coincident_lower_d_calc;
+
   /// Quadrature point index
   unsigned int _qp;
 

framework/include/base/Assembly.h

@@ -1906,6 +1906,11 @@ class Assembly
    */
   void havePRefinement(const std::vector<FEFamily> & disable_p_refinement_for_families);
 
+  /**
+   * Set the current lower dimensional element. This may be null
+   */
+  void setCurrentLowerDElem(const Elem * const lower_d_elem);
+
 private:
   /**
    * Just an internal helper function to reinit the volume FE objects.
@@ -3153,3 +3158,9 @@ Assembly::assignDisplacements(
 {
   _disp_numbers_and_directions = std::move(disp_numbers_and_directions);
 }
+
+inline void
+Assembly::setCurrentLowerDElem(const Elem * const lower_d_elem)
+{
+  _current_lower_d_elem = lower_d_elem;
+}

framework/include/problems/FEProblemBase.h

@@ -2144,6 +2144,8 @@ class FEProblemBase : public SubProblem, public Restartable
    */
   bool identifyVariableGroupsInNL() const { return _identify_variable_groups_in_nl; }
 
+  virtual void setCurrentLowerDElem(const Elem * const lower_d_elem, const THREAD_ID tid) override;
+
 protected:
   /// Create extra tagged vectors and matrices
   void createTagVectors();

framework/include/problems/SubProblem.h

@@ -939,6 +939,11 @@ class SubProblem : public Problem
    */
   [[nodiscard]] bool havePRefinement() const { return _have_p_refinement; }
 
+  /**
+   * Set the current lower dimensional element. This can be null
+   */
+  virtual void setCurrentLowerDElem(const Elem * const lower_d_elem, const THREAD_ID tid);
+
 protected:
   /**
    * Helper function called by getVariable that handles the logic for

framework/include/variables/MooseVariableBase.h

@@ -178,6 +178,11 @@ class MooseVariableBase : public MooseObject,
    */
   bool isArray() const { return _is_array; }
 
+  /**
+   * @return whether this variable lives on lower dimensional blocks
+   */
+  bool isLowerD() const { return _is_lower_d; }
+
 protected:
   /**
    * @returns whether we should insert derivatives
@@ -237,6 +242,9 @@ class MooseVariableBase : public MooseObject,
 
   /// Whether this is an array variable
   const bool _is_array;
+
+  /// Whether this variable lives on lower dimensional blocks
+  const bool _is_lower_d;
 };
 
 inline void

framework/include/variables/MooseVariableFE.h

@@ -511,7 +511,8 @@ class MooseVariableFE : public MooseVariableField<OutputType>
   /**
    * Set local DOF values and evaluate the values on quadrature points
    */
-  void setDofValues(const DenseVector<OutputData> & values) override;
+  virtual void setDofValues(const DenseVector<OutputData> & values) override;
+  virtual void setLowerDofValues(const DenseVector<OutputData> & values) override;
 
   /**
    * Write a nodal value to the passed-in solution vector
@@ -551,18 +552,23 @@ class MooseVariableFE : public MooseVariableField<OutputType>
    * @return Variable value
    */
   OutputData getElementalValueOlder(const Elem * elem, unsigned int idx = 0) const;
+
   /**
    * Set the current local DOF values to the input vector
    */
-  void insert(NumericVector<Number> & residual) override;
+  virtual void insert(NumericVector<Number> & vector) override;
+  virtual void insertLower(NumericVector<Number> & vector) override;
+
   /**
    * Add the current local DOF values to the input vector
    */
-  void add(NumericVector<Number> & residual) override;
+  virtual void add(NumericVector<Number> & vector) override;
+
   /**
    * Add passed in local DOF values onto the current solution
    */
   void addSolution(const DenseVector<Number> & v);
+
   /**
    * Add passed in local neighbor DOF values onto the current solution
    */

framework/include/variables/MooseVariableFV.h

@@ -379,7 +379,8 @@ class MooseVariableFV : public MooseVariableField<OutputType>
   /**
    * Set local DOF values and evaluate the values on quadrature points
    */
-  void setDofValues(const DenseVector<OutputData> & values) override;
+  virtual void setDofValues(const DenseVector<OutputData> & values) override;
+  virtual void setLowerDofValues(const DenseVector<OutputData> & values) override;
 
   /// Get the current value of this variable on an element
   /// @param[in] elem   Element at which to get value
@@ -397,8 +398,9 @@ class MooseVariableFV : public MooseVariableField<OutputType>
   /// @return Variable value
   OutputData getElementalValueOlder(const Elem * elem, unsigned int idx = 0) const;
 
-  virtual void insert(NumericVector<Number> & residual) override;
-  virtual void add(NumericVector<Number> & residual) override;
+  virtual void insert(NumericVector<Number> & vector) override;
+  virtual void insertLower(NumericVector<Number> & vector) override;
+  virtual void add(NumericVector<Number> & vector) override;
 
   const DoFValue & dofValues() const override;
   const DoFValue & dofValuesOld() const override;
@@ -684,6 +686,11 @@ class MooseVariableFV : public MooseVariableField<OutputType>
   /// in \p getDirichletBC
   std::unordered_map<BoundaryID, const FVDirichletBCBase *> _boundary_id_to_dirichlet_bc;
 
+  /**
+   * Emit an error message for unsupported lower-d ops
+   */
+  [[noreturn]] void lowerDError() const;
+
 protected:
   /// A cache for storing gradients on elements
   mutable std::unordered_map<const Elem *, VectorValue<ADReal>> _elem_to_grad;
@@ -819,7 +826,14 @@ template <typename OutputType>
 const typename MooseVariableFV<OutputType>::FieldVariablePhiValue &
 MooseVariableFV<OutputType>::phiLower() const
 {
-  mooseError("Not defined for finite volume variables");
+  lowerDError();
+}
+
+template <typename OutputType>
+void
+MooseVariableFV<OutputType>::lowerDError() const
+{
+  mooseError("Lower dimensional element support not implemented for finite volume variables");
 }
 
 template <>

framework/include/variables/MooseVariableField.h

@@ -304,6 +304,12 @@ class MooseVariableField : public MooseVariableFieldBase,
    */
   virtual void setDofValues(const DenseVector<OutputData> & values) = 0;
 
+  /**
+   * Set local DOF values for a lower dimensional element and evaluate the values on quadrature
+   * points
+   */
+  virtual void setLowerDofValues(const DenseVector<OutputData> & values) = 0;
+
   /**
    * Whether or not this variable is actually using the shape function value.
    *

framework/include/variables/MooseVariableFieldBase.h

@@ -173,8 +173,21 @@ class MooseVariableFieldBase : public MooseVariableBase
 
   virtual unsigned int numberOfDofsNeighbor() = 0;
 
-  virtual void insert(NumericVector<Number> & residual) = 0;
-  virtual void add(NumericVector<Number> & residual) = 0;
+  /**
+   * Insert the currently cached degree of freedom values into the provided \p vector
+   */
+  virtual void insert(NumericVector<Number> & vector) = 0;
+
+  /**
+   * Insert the currently cached degree of freedom values for a lower-dimensional element into the
+   * provided \p vector
+   */
+  virtual void insertLower(NumericVector<Number> & vector) = 0;
+
+  /**
+   * Add the currently cached degree of freedom values into the provided \p vector
+   */
+  virtual void add(NumericVector<Number> & vector) = 0;
 
   /**
    * Return phi size

framework/src/auxkernels/AuxKernel.C

@@ -138,12 +138,15 @@ AuxKernelTempl<ComputeValueType>::AuxKernelTempl(const InputParameters & paramet
 
     _current_node(_assembly.node()),
     _current_boundary_id(_assembly.currentBoundaryID()),
-    _solution(_aux_sys.solution())
+    _solution(_aux_sys.solution()),
+
+    _current_lower_d_elem(_assembly.lowerDElem()),
+    _coincident_lower_d_calc(_bnd && !isNodal() && _var.isLowerD())
 {
   addMooseVariableDependency(&_var);
   _supplied_vars.insert(parameters.get<AuxVariableName>("variable"));
 
-  if (_bnd && !isNodal() && _check_boundary_restricted)
+  if (_bnd && !isNodal() && !_coincident_lower_d_calc && _check_boundary_restricted)
   {
     // when the variable is elemental and this aux kernel operates on boundaries,
     // we need to check that no elements are visited more than once through visiting
@@ -274,6 +277,16 @@ AuxKernelTempl<RealVectorValue>::setDofValueHelper(const RealVectorValue &)
   mooseError("Not implemented");
 }
 
+template <typename ComputeValueType>
+void
+AuxKernelTempl<ComputeValueType>::insert()
+{
+  if (_coincident_lower_d_calc)
+    _var.insertLower(_aux_sys.solution());
+  else
+    _var.insert(_aux_sys.solution());
+}
+
 template <typename ComputeValueType>
 void
 AuxKernelTempl<ComputeValueType>::compute()
@@ -282,6 +295,9 @@ AuxKernelTempl<ComputeValueType>::compute()
 
   if (isNodal()) /* nodal variables */
   {
+    mooseAssert(!_coincident_lower_d_calc,
+                "Nodal evaluations are point evaluations. We don't have to concern ourselves with "
+                "coincidence of lower-d blocks and higher-d faces because they share nodes");
     if (_var.isNodalDefined())
     {
       _qp = 0;
@@ -292,7 +308,18 @@ AuxKernelTempl<ComputeValueType>::compute()
   }
   else /* elemental variables */
   {
-    _n_local_dofs = _var.numberOfDofs();
+    _n_local_dofs = _coincident_lower_d_calc ? _var.dofIndicesLower().size() : _var.numberOfDofs();
+
+    if (_coincident_lower_d_calc)
+    {
+      static const std::string lower_error = "Make sure that the lower-d variable lives on a "
+                                             "lower-d block that is a superset of the boundary";
+      if (!_current_lower_d_elem)
+        mooseError("No lower-dimensional element. ", lower_error);
+      if (!_n_local_dofs)
+        mooseError("No degrees of freedom. ", lower_error);
+    }
+
     if (_n_local_dofs == 1) /* p0 */
     {
       ComputeValueType value = 0;
@@ -302,10 +329,22 @@ AuxKernelTempl<ComputeValueType>::compute()
       if (_var.isFV())
         setDofValueHelper(value);
       else
+      {
         // update the variable data referenced by other kernels.
         // Note that this will update the values at the quadrature points too
         // (because this is an Elemental variable)
-        _var.setNodalValue(value);
+        if (_coincident_lower_d_calc)
+        {
+          _local_sol.resize(1);
+          if constexpr (std::is_same<Real, ComputeValueType>::value)
+            _local_sol(0) = value;
+          else
+            mooseAssert(false, "We should not enter the single dof branch with a vector variable");
+          _var.setLowerDofValues(_local_sol);
+        }
+        else
+          _var.setNodalValue(value);
+      }
     }
     else /* high-order */
     {
@@ -314,14 +353,16 @@ AuxKernelTempl<ComputeValueType>::compute()
       _local_ke.resize(_n_local_dofs, _n_local_dofs);
       _local_ke.zero();
 
+      const auto & test = _coincident_lower_d_calc ? _var.phiLower() : _test;
+
       // assemble the local mass matrix and the load
-      for (unsigned int i = 0; i < _test.size(); i++)
+      for (unsigned int i = 0; i < test.size(); i++)
         for (_qp = 0; _qp < _qrule->n_points(); _qp++)
         {
-          ComputeValueType t = _JxW[_qp] * _coord[_qp] * _test[i][_qp];
+          ComputeValueType t = _JxW[_qp] * _coord[_qp] * test[i][_qp];
           _local_re(i) += t * computeValue();
-          for (unsigned int j = 0; j < _test.size(); j++)
-            _local_ke(i, j) += t * _test[j][_qp];
+          for (unsigned int j = 0; j < test.size(); j++)
+            _local_ke(i, j) += t * test[j][_qp];
         }
       // mass matrix is always SPD but in case of boundary restricted, it will be rank deficient
       _local_sol.resize(_n_local_dofs);
@@ -330,7 +371,7 @@ AuxKernelTempl<ComputeValueType>::compute()
       else
         _local_ke.cholesky_solve(_local_re, _local_sol);
 
-      _var.setDofValues(_local_sol);
+      _coincident_lower_d_calc ? _var.setLowerDofValues(_local_sol) : _var.setDofValues(_local_sol);
     }
   }
 }

framework/src/loops/ComputeElemAuxBcsThread.C

@@ -82,6 +82,7 @@ ComputeElemAuxBcsThread<AuxKernelType>::operator()(const ConstBndElemRange & ran
         _fe_problem.reinitElemFace(elem, side, boundary_id, _tid);
 
         const Elem * lower_d_elem = _fe_problem.mesh().getLowerDElem(elem, side);
+        _fe_problem.setCurrentLowerDElem(lower_d_elem, _tid);
         if (lower_d_elem)
           _fe_problem.reinitLowerDElem(lower_d_elem, _tid);
 
@@ -126,7 +127,7 @@ ComputeElemAuxBcsThread<AuxKernelType>::operator()(const ConstBndElemRange & ran
         for (const auto & aux : iter->second)
         {
           aux->compute();
-          aux->variable().insert(_aux_sys.solution());
+          aux->insert();
         }
 
         if (_need_materials)

framework/src/meshgenerators/RefineSidesetGenerator.C

@@ -52,6 +52,9 @@ RefineSidesetGenerator::RefineSidesetGenerator(const InputParameters & parameter
   if (_boundaries.size() != _refinement.size())
     paramError("refinement",
                "The boundaries and refinement parameter vectors should be the same size");
+  if (_boundaries.size() != _boundary_side.size())
+    paramError("boundary_side",
+               "The boundaries and boundary_side parameter vectors should be the same size");
 }
 
 std::unique_ptr<MeshBase>