Add weight functions

framework/doc/content/source/userobject/RadialAverage.md

@@ -6,9 +6,10 @@ Given a material property and a radius for averaging, the RadialAverage object
 computes the spatial average value of the property over the radius.
 
 ## Applications
+
 This can be used for nonlocal damage models in the `TensorMechanics` module
 where the damage_index that is used for computing the damage stress is average
-over a certain radius. This can help alleviate mesh sensitivity in certain
+over a certain radius $r_0$. This can help alleviate mesh sensitivity in certain
 cases. This can be accomplished by running the RadialAverage object on a local
 damage material property. Then using the `NonlocalDamage` model in conjunction
 with the `ComputeDamageStress` the damage index used for updating the stress is
@@ -20,7 +21,7 @@ The RadialAverage user object is derived from `ElementUserObject` and
 works in two stages.
 
 1. In the element loop (in the `execute()` method) a list of all quadrature
-   points, their locations, indices, and selected variable value is compiled.
+   points, their locations, indices, and selected material property value is compiled.
 
 2. In the `finalize()` method
 
@@ -36,6 +37,11 @@ works in two stages.
     4. the results from the search are used to spatially averaged
 
 
+The [!param](/UserObjects/RadialAverage/weights) parameter determines the distance
+based weight function to be used in the averaging process. `constant` assigns an equal weight
+to all material points, `linear` weights each point by $r_0-r$ (a linear fall-off with distance),
+and `cosine` weights each point by $1+\cos (\frac r{r_0}\pi}$.
+
 !syntax parameters /UserObjects/RadialAverage
 
 !syntax inputs /UserObjects/RadialAverage

framework/include/userobject/RadialAverage.h

@@ -72,10 +72,16 @@ class RadialAverage : public ElementUserObject
 protected:
   void updateCommunicationLists();
 
-  /// material name to get averaged
-  std::string _averaged_material_name;
+  /// distance based weight function
+  enum class WeightsType
+  {
+    CONSTANT,
+    LINEAR,
+    COSINE
+  } _weights_type;
+
   /// material to be averaged
-  const MaterialProperty<Real> & _averaged_material;
+  const MaterialProperty<Real> & _prop;
 
   /// cut-off radius
   const Real _radius;
@@ -111,6 +117,9 @@ class RadialAverage : public ElementUserObject
   std::vector<std::set<std::size_t>> _communication_lists;
   bool _update_communication_lists;
 
+  /// processors to send (potentially empty) data to
+  std::vector<processor_id_type> _candidate_procs;
+
   processor_id_type _my_pid;
 
   //@{ PerfGraph identifiers

framework/src/userobject/RadialAverage.C

@@ -34,8 +34,11 @@ InputParameters
 RadialAverage::validParams()
 {
   InputParameters params = ElementUserObject::validParams();
-  params.addClassDescription("Perform a radial equal weight average of a material property");
-  params.addRequiredParam<std::string>("material_name", "Name of the material to average.");
+  params.addClassDescription("Perform a radial average of a material property");
+  params.addRequiredParam<MaterialPropertyName>("prop_name",
+                                                "Name of the material property to average");
+  MooseEnum weights_type("constant linear cosine", "linear");
+  params.addRequiredParam<MooseEnum>("weights", weights_type, "Distance based weight function");
   params.addRequiredParam<Real>("radius", "Cut-off radius for the averaging");
   params.addRangeCheckedParam<Real>(
       "padding",
@@ -57,8 +60,8 @@ RadialAverage::validParams()
 
 RadialAverage::RadialAverage(const InputParameters & parameters)
   : ElementUserObject(parameters),
-    _averaged_material_name(getParam<std::string>("material_name")),
-    _averaged_material(getMaterialProperty<Real>(_averaged_material_name)),
+    _weights_type(getParam<MooseEnum>("weights").getEnum<WeightsType>()),
+    _prop(getMaterialProperty<Real>("prop_name")),
     _radius(getParam<Real>("radius")),
     _padding(getParam<Real>("padding")),
     _update_communication_lists(false),
@@ -89,7 +92,7 @@ RadialAverage::execute()
 
   // collect all QP data
   for (const auto qp : make_range(_qrule->n_points()))
-    _qp_data.emplace_back(_q_point[qp], id, qp, _JxW[qp] * _coord[qp], _averaged_material[qp]);
+    _qp_data.emplace_back(_q_point[qp], id, qp, _JxW[qp] * _coord[qp], _prop[qp]);
 
   // make sure the result map entry for the current element is sized correctly
   auto i = _average.find(id);
@@ -120,24 +123,18 @@ RadialAverage::finalize()
         libMesh::n_threads() > 1)
       updateCommunicationLists();
 
-    // sparse send data (processor ID,)
-    std::vector<std::size_t> non_zero_comm;
-    for (auto i = beginIndex(_communication_lists); i < _communication_lists.size(); ++i)
-      if (!_communication_lists[i].empty())
-        non_zero_comm.push_back(i);
-
     // data structures for sparse point to point communication
-    std::vector<std::vector<QPData>> send(non_zero_comm.size());
-    std::vector<Parallel::Request> send_requests(non_zero_comm.size());
+    std::vector<std::vector<QPData>> send(_candidate_procs.size());
+    std::vector<Parallel::Request> send_requests(_candidate_procs.size());
     Parallel::MessageTag send_tag = _communicator.get_unique_tag(4711);
     std::vector<QPData> receive;
 
     const auto item_type = TIMPI::StandardType<QPData>(&(_qp_data[0]));
 
     // fill buffer and send structures
-    for (auto i = beginIndex(non_zero_comm); i < non_zero_comm.size(); ++i)
+    for (const auto i : index_range(_candidate_procs))
     {
-      const auto pid = non_zero_comm[i];
+      const auto pid = _candidate_procs[i];
       const auto & list = _communication_lists[pid];
 
       // fill send buffer for transfer to pid
@@ -150,8 +147,12 @@ RadialAverage::finalize()
     }
 
     // receive messages - we assume that we receive as many messages as we send!
-    for (auto i = beginIndex(non_zero_comm); i < non_zero_comm.size(); ++i)
+    // bounding box overlapp is transitive, but data exhange between overlapping procs could still
+    // be unidirectional!
+    for (const auto i : index_range(_candidate_procs))
     {
+      libmesh_ignore(i);
+
       // inspect incoming message
       Parallel::Status status(_communicator.probe(Parallel::any_source, send_tag));
       const auto source_pid = TIMPI::cast_int<processor_id_type>(status.source());
@@ -276,14 +277,14 @@ RadialAverage::updateCommunicationLists()
     bbs.emplace_back(pp.first - rpoint, pp.second + rpoint);
 
   // get candidate processors (overlapping bounding boxes)
-  std::vector<processor_id_type> candidate_procs;
+  _candidate_procs.clear();
   for (const auto pid : index_range(bbs))
     if (pid != _my_pid && bbs[pid].intersects(mypp))
-      candidate_procs.push_back(pid);
+      _candidate_procs.push_back(pid);
 
   // go over all boundary data items and send them to the proc they overlap with
   for (const auto i : _boundary_data_indices)
-    for (const auto pid : candidate_procs)
+    for (const auto pid : _candidate_procs)
       if (bbs[pid].contains_point(_qp_data[i]._q_point))
         _communication_lists[pid].insert(i);
 

framework/src/userobject/ThreadedRadialAverageLoop.C

@@ -26,6 +26,7 @@ ThreadedRadialAverageLoop::operator()(const QPDataRange & qpdata_range)
   const auto radius = _radavg._radius;
   const auto & qp_data = _radavg._qp_data;
   const auto & kd_tree = _radavg._kd_tree;
+  const auto & weights_type = _radavg._weights_type;
 
   // tree search data structures
   std::vector<std::pair<std::size_t, Real>> ret_matches;
@@ -52,11 +53,26 @@ ThreadedRadialAverageLoop::operator()(const QPDataRange & qpdata_range)
     std::size_t n_result =
         kd_tree->radiusSearch(&(local_qp._q_point(0)), radius * radius, ret_matches, search_params);
     Real total_vol = 0.0;
+    Real weight = 1.0;
     for (std::size_t j = 0; j < n_result; ++j)
     {
       const auto & other_qp = qp_data[ret_matches[j].first];
-      sum += other_qp._value * other_qp._volume;
-      total_vol += other_qp._volume;
+      switch (weights_type)
+      {
+        case RadialAverage::WeightsType::CONSTANT:
+          break;
+
+        case RadialAverage::WeightsType::LINEAR:
+          weight = radius - std::sqrt(ret_matches[j].second);
+          break;
+
+        case RadialAverage::WeightsType::COSINE:
+          weight = std::cos(std::sqrt(ret_matches[j].second) / radius * libMesh::pi) + 1.0;
+          break;
+      }
+
+      sum += other_qp._value * other_qp._volume * weight;
+      total_vol += other_qp._volume * weight;
     }
     sum /= total_vol;
   }

modules/tensor_mechanics/test/tests/scalar_material_damage/ad_nonlocal_scalar_damage.i

@@ -60,9 +60,9 @@
 [UserObjects]
   [ele_avg]
     type = RadialAverage
-    material_name = local_damage_reg
+    prop_name = local_damage_reg
+    weights = constant
     execute_on = "INITIAL timestep_end"
-    block = 0
     radius = 0.55
   []
 []

modules/tensor_mechanics/test/tests/scalar_material_damage/nonlocal_scalar_damage.i

@@ -59,9 +59,9 @@
 [UserObjects]
   [ele_avg]
     type = RadialAverage
-    material_name = local_damage
+    prop_name = local_damage
+    weights = constant
     execute_on = "INITIAL timestep_end"
-    block = 0
     radius = 0.55
   []
 []

test/tests/userobjects/radial_average/time_changing_test.i → test/tests/userobjects/radial_average/test.i

@@ -18,13 +18,14 @@
 [AuxVariables]
   [non_local_material]
     family = MONOMIAL
+    order = SECOND
   []
 []
 
 [AuxKernels]
   [non_local]
     type = RadialAverageAux
-    average_UO = ele_avg
+    average_UO = average
     variable = non_local_material
     execute_on = 'INITIAL TIMESTEP_BEGIN'
   []
@@ -46,11 +47,11 @@
 []
 
 [UserObjects]
-  [ele_avg]
+  [average]
     type = RadialAverage
-    material_name = local
+    prop_name = local
+    weights = constant
     execute_on = "INITIAL timestep_end"
-    block = 0
     radius = 0.3
   []
 []

test/tests/userobjects/radial_average/tests

@@ -1,13 +1,36 @@
 [Tests]
-  issues = '#21786'
+  issues = '#21786 #21944 #21963'
   design = '/RadialAverage.md'
 
-  [radial_average]
-    type = 'Exodiff'
-    input = 'time_changing_test.i'
-    exodiff = 'time_changing_test_out.e'
-    recover = false
-    requirement = 'The system shall will compute a radial average of a material that changes over time.'
-    abs_zero = 1e-10
+  [weight]
+    requirement = 'The system shall will compute a radial average of a material property that changes over time with '
+    [constant]
+      type = 'Exodiff'
+      input = 'test.i'
+      cli_args = 'Outputs/file_base=constant UserObjects/average/weights=constant'
+      exodiff = 'constant.e'
+      recover = false
+      abs_zero = 1e-10
+      detail = 'an equal weight for all material points'
+    []
+    [linear]
+      type = 'Exodiff'
+      input = 'test.i'
+      cli_args = 'Outputs/file_base=linear UserObjects/average/weights=linear'
+      exodiff = 'linear.e'
+      recover = false
+      abs_zero = 1e-10
+      detail = 'a weight factor that falls off linearly to zero with distance'
+    []
+    [cosine]
+      type = 'Exodiff'
+      input = 'test.i'
+      cli_args = 'Outputs/file_base=cosine UserObjects/average/weights=cosine'
+      exodiff = 'cosine.e'
+      recover = false
+      abs_zero = 1e-10
+      detail = 'a weight factor that falls off with a cosine profile with distance'
+    []
   []
 []
+