Add side Qp support to SideIntegralFunctorPP

refs idaholab#16099

framework/include/postprocessors/SideIntegralFunctorPostprocessor.h

@@ -21,7 +21,7 @@ typedef SideIntegralFunctorPostprocessorTempl<true> ADSideIntegralFunctorPostpro
  * This postprocessor computes a surface integral of the specified functor
  *
  * Note that specializations of this integral are possible by deriving from this
- * class and overriding computeQpIntegral().
+ * class and overriding computeQpIntegral() or computeFaceInfoIntegral()
  */
 template <bool is_ad>
 class SideIntegralFunctorPostprocessorTempl : public SideIntegralPostprocessor
@@ -39,11 +39,13 @@ class SideIntegralFunctorPostprocessorTempl : public SideIntegralPostprocessor
    */
   virtual Real computeFaceInfoIntegral(const FaceInfo * fi) override;
 
-  /// Quadrature point integration of functors is not implemented
-  Real computeQpIntegral() override
-  {
-    mooseError("Functor side integrals are implemented over FaceInfos not quadrature points");
-  };
+  Real computeQpIntegral() override;
+
+  /// Check if the functor and the prefactor are defined on the block right by the subdomain
+  bool checkFunctorDefinedOnSideBlock() const;
+
+  /// Error with a helpful message if the functor is not defined on the block by the sideset
+  void errorFunctorNotDefinedOnSideBlock() const;
 
   /// Functor being integrated
   const Moose::Functor<GenericReal<is_ad>> & _functor;

framework/src/postprocessors/SideIntegralFunctorPostprocessor.C

@@ -28,6 +28,8 @@ SideIntegralFunctorPostprocessorTempl<is_ad>::validParams()
                         "the sideset, allows to skip the parts where it is not defined. Please "
                         "keep in mind that if the sideset is defined with the wrong normal, this "
                         "may allow to skip the entire integral.");
+  MooseEnum functor_args("qp face", "face");
+  params.addParam<MooseEnum>("functor_argument", functor_args, "Location to evaluated functors at");
   params.addClassDescription(
       "Computes a surface integral of the specified functor, using the "
       "single-sided face argument, which usually means that the functor will be"
@@ -44,8 +46,7 @@ SideIntegralFunctorPostprocessorTempl<is_ad>::SideIntegralFunctorPostprocessorTe
     _prefactor(getFunctor<GenericReal<is_ad>>("prefactor")),
     _partial_integral(getParam<bool>("restrict_to_functors_domain"))
 {
-  // Only support face info integration for now
-  _qp_integration = false;
+  _qp_integration = (getParam<MooseEnum>("functor_argument") == "qp");
 }
 
 template <bool is_ad>
@@ -57,43 +58,75 @@ SideIntegralFunctorPostprocessorTempl<is_ad>::computeFaceInfoIntegral(const Face
   // We wont allow that case, unless explicitly requested by the user,
   // as this means the sideset is reversed, but we will allow the case where
   // both sides are defined
-  bool has_elem;
+  const bool has_elem = checkFunctorDefinedOnSideBlock();
+
+  if (has_elem)
+  {
+    Moose::SingleSidedFaceArg ssf = {
+        fi, Moose::FV::LimiterType::CentralDifference, true, false, _current_elem->subdomain_id()};
+    return MetaPhysicL::raw_value(_prefactor(ssf) * _functor(ssf));
+  }
+  else
+  {
+    if (!_partial_integral)
+      errorFunctorNotDefinedOnSideBlock();
+    return 0;
+  }
+}
+
+template <bool is_ad>
+Real
+SideIntegralFunctorPostprocessorTempl<is_ad>::computeQpIntegral()
+{
+  const bool has_elem = checkFunctorDefinedOnSideBlock();
+
+  if (has_elem)
+  {
+    Moose::ElemSideQpArg elem_side_qp = {_current_elem, _current_side, _qp, _qrule};
+    return MetaPhysicL::raw_value(_prefactor(elem_side_qp) * _functor(elem_side_qp));
+  }
+  else
+  {
+    if (!_partial_integral)
+      errorFunctorNotDefinedOnSideBlock();
+    return 0;
+  }
+}
+
+template <bool is_ad>
+bool
+SideIntegralFunctorPostprocessorTempl<is_ad>::checkFunctorDefinedOnSideBlock() const
+{
   if (_functor.hasBlocks(_current_elem->subdomain_id()) &&
       _prefactor.hasBlocks(_current_elem->subdomain_id()))
-    has_elem = true;
+    return true;
   else
   {
 #ifndef NDEBUG
-    if (fi->neighborPtr())
+    const auto neighbor_ptr = _current_elem->neighbor_ptr(_current_side);
+    if (neighbor_ptr)
     {
       mooseAssert(_functor.hasBlocks(_current_elem->subdomain_id()) ||
-                      _functor.hasBlocks(fi->neighborPtr()->subdomain_id()),
+                      _functor.hasBlocks(neighbor_ptr->subdomain_id()),
                   "Functor should be defined on at least one side of the boundary");
       mooseAssert(_prefactor.hasBlocks(_current_elem->subdomain_id()) ||
-                      _prefactor.hasBlocks(fi->neighborPtr()->subdomain_id()),
+                      _prefactor.hasBlocks(neighbor_ptr->subdomain_id()),
                   "Prefactor should be defined on at least one side of the boundary");
     }
 #endif
-    has_elem = false;
+    return false;
   }
+}
 
-  if (has_elem)
-  {
-    Moose::SingleSidedFaceArg ssf = {
-        fi, Moose::FV::LimiterType::CentralDifference, true, false, _current_elem->subdomain_id()};
-    return MetaPhysicL::raw_value(_prefactor(ssf) * _functor(ssf));
-  }
-  else
-  {
-    if (!_partial_integral)
-      paramError("boundary",
-                 "Functor " + _functor.functorName() + " (or prefactor " +
-                     _prefactor.functorName() + ") is not defined on block " +
-                     std::to_string(_current_elem->subdomain_id()) +
-                     ". Is the functor defined along the whole sideset? "
-                     "Are the sidesets in 'boundary' all oriented correctly?");
-    return 0;
-  }
+template <bool is_ad>
+void
+SideIntegralFunctorPostprocessorTempl<is_ad>::errorFunctorNotDefinedOnSideBlock() const
+{
+  paramError("boundary",
+             "Functor " + _functor.functorName() + " (or prefactor " + _prefactor.functorName() +
+                 ") is not defined on block " + std::to_string(_current_elem->subdomain_id()) +
+                 ". Is the functor defined along the whole sideset? "
+                 "Are the sidesets in 'boundary' all oriented correctly?");
 }
 
 template class SideIntegralFunctorPostprocessorTempl<false>;