Add point reduction algorithm (idaholab#14220)

framework/doc/content/source/functions/CoarsendPiecewiseLinear.md

@@ -1,10 +1,10 @@
 # PiecewiseLinear
 
-!syntax description /Functions/CoarsenedPiecewiseLinear
+!syntax description /Functions/CoarsendPiecewiseLinear
 
 ## Description
 
-The `CoarsenedPiecewiseLinear` performs preprocessing and linear interpolation
+The `CoarsendPiecewiseLinear` performs preprocessing and linear interpolation
 on an x/y data set. The object  acts like
 [`PiecewiseLinear`](/PiecewiseLinear.md)  except that it reduces the number of
 function point at the start of the simulation. It uses the
@@ -15,8 +15,8 @@ for data reduction.
 
 !listing test/tests/misc/check_error/function_file_test1.i block=Functions
 
-!syntax parameters /Functions/CoarsenedPiecewiseLinear
+!syntax parameters /Functions/CoarsendPiecewiseLinear
 
-!syntax inputs /Functions/CoarsenedPiecewiseLinear
+!syntax inputs /Functions/CoarsendPiecewiseLinear
 
-!syntax children /Functions/CoarsenedPiecewiseLinear
+!syntax children /Functions/CoarsendPiecewiseLinear

framework/doc/content/source/vectorpostprocessors/PiecewiseFunctionTabulate.md

@@ -0,0 +1,14 @@
+# PiecewiseFunctionTabulate
+
+!syntax description /VectorPostprocessors/PiecewiseFunctionTabulate
+
+This object can be used to check the function generated by a
+[`CoarsenedPicewiseLinear`](/CoarsenedPicewiseLinear.md) function.
+
+!syntax parameters /VectorPostprocessors/PiecewiseFunctionTabulate
+
+!syntax inputs /VectorPostprocessors/PiecewiseFunctionTabulate
+
+!syntax children /VectorPostprocessors/PiecewiseFunctionTabulate
+
+!bibtex bibliography

framework/include/utils/PointReduction.h

@@ -0,0 +1,29 @@
+//* 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 "Moose.h"
+#include <vector>
+
+namespace PointReduction
+{
+
+typedef std::pair<Real, Real> FunctionNode;
+typedef std::vector<FunctionNode> FunctionNodeList;
+
+/// compute the perpendicular distance of a point from a line defined by begin and end points
+Real perpendicularDistance(const FunctionNode & point,
+                           const FunctionNode & begin,
+                           const FunctionNode & end);
+
+/// return a pruned function node list using the Ramer-Douglas-Peucker algorithm
+FunctionNodeList douglasPeucker(const FunctionNodeList &, Real epsilon);
+
+} // namespace PointReduction

framework/include/vectorpostprocessors/PiecewiseFunctionTabulate.h

@@ -0,0 +1,37 @@
+//* 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 "GeneralVectorPostprocessor.h"
+
+class PiecewiseFunctionTabulate;
+class PiecewiseBase;
+
+template <>
+InputParameters validParams<PiecewiseFunctionTabulate>();
+
+/**
+ * Tabulate the function nodes of a piecewise function, such as PiecewiseLinear or
+ * PiecewiseConstant
+ */
+class PiecewiseFunctionTabulate : public GeneralVectorPostprocessor
+{
+public:
+  PiecewiseFunctionTabulate(const InputParameters & parameters);
+
+  virtual void initialize() override {}
+  virtual void execute() override;
+  virtual void finalize() override {}
+
+protected:
+  const PiecewiseBase * _piecewise_function;
+  VectorPostprocessorValue & _x_col;
+  VectorPostprocessorValue & _y_col;
+};

framework/src/functions/CoarsendPiecewiseLinear.C

@@ -8,6 +8,8 @@
 //* https://www.gnu.org/licenses/lgpl-2.1.html
 
 #include "CoarsendPiecewiseLinear.h"
+#include "PointReduction.h"
+#include "TimedPrint.h"
 
 registerMooseObject("MooseApp", CoarsendPiecewiseLinear);
 
@@ -18,11 +20,48 @@ validParams<CoarsendPiecewiseLinear>()
   InputParameters params = validParams<PiecewiseLinearBase>();
   params.addClassDescription("Perform a point reduction of the tabulated data upon initialization, "
                              "then evaluate using a linear interpolation.");
+  params.addRequiredParam<Real>(
+      "epsilon",
+      "Significant distance in the function below which points are considered removable noise");
+  params.addParam<Real>("y_scale",
+                        1.0,
+                        "Scaling factor to apply to the function nodes for the purpose of "
+                        "computing distances in the Douglas-Peucker point reduction algorithm. "
+                        "This permits shifting the weight between x and y-direction distances.");
+  params.addParam<Real>("x_scale",
+                        1.0,
+                        "Scaling factor to apply to the function nodes for the purpose of "
+                        "computing distances in the Douglas-Peucker point reduction algorithm. "
+                        "This permits shifting the weight between x and y-direction distances.");
   return params;
 }
 
 CoarsendPiecewiseLinear::CoarsendPiecewiseLinear(const InputParameters & parameters)
   : PiecewiseLinearBase(parameters)
 {
+  const Real x_scale = getParam<Real>("x_scale");
+  const Real y_scale = getParam<Real>("y_scale");
+  const Real epsilon = getParam<Real>("epsilon");
+
+  // create vector of pairs
+  PointReduction::FunctionNodeList list;
+  list.reserve(_raw_x.size());
+  for (MooseIndex(_raw_x) i = 0; i < _raw_x.size(); ++i)
+    list.emplace_back(_raw_x[i] * x_scale, _raw_y[i] * y_scale);
+
+  // point reduction
+  _console << "Reduced size for function '" << name() << "' from " << list.size();
+  list = PointReduction::douglasPeucker(list, epsilon);
+  _console << " to " << list.size() << " points.\n";
+
+  // unpack vector of pairs
+  _raw_x.resize(list.size());
+  _raw_y.resize(list.size());
+  for (MooseIndex(list) i = 0; i < list.size(); ++i)
+  {
+    _raw_x[i] = list[i].first / x_scale;
+    _raw_y[i] = list[i].second / y_scale;
+  }
+
   buildInterpolation();
 }

framework/src/utils/PointReduction.C

@@ -0,0 +1,98 @@
+//* 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 "PointReduction.h"
+#include "MooseError.h"
+
+#include <algorithm>
+#include <cmath>
+
+namespace PointReduction
+{
+
+Real
+sqr(Real a)
+{
+  return a * a;
+}
+
+Real
+perpendicularDistance(const FunctionNode & node,
+                      const FunctionNode & begin,
+                      const FunctionNode & end)
+{
+  const Real x0 = node.first;
+  const Real y0 = node.second;
+  const Real x1 = begin.first;
+  const Real y1 = begin.second;
+  const Real x2 = end.first;
+  const Real y2 = end.second;
+
+  const Real denom = std::sqrt(sqr(y2 - y1) + sqr(x2 - x1));
+  mooseAssert(denom > 0, "Line begin and end points bust not be the same");
+
+  return std::abs((y2 - y1) * x0 - (x2 - x1) * y0 + x2 * y1 - y2 * x1) / denom;
+}
+
+void
+douglasPeuckerRecurse(const FunctionNodeList & list,
+                      const Real epsilon,
+                      std::vector<bool> & keep,
+                      std::size_t begin,
+                      std::size_t end)
+{
+  // Find the point with the maximum distance
+  Real dmax = 0.0;
+  std::size_t index = 0;
+
+  for (std::size_t i = begin; i <= end; ++i)
+    if (keep[i])
+    {
+      const Real d = perpendicularDistance(list[i], list[begin], list[end]);
+      if (d > dmax)
+      {
+        index = i;
+        dmax = d;
+      }
+    }
+
+  // If max distance is greater than epsilon, recursively simplify
+  if (dmax > epsilon)
+  {
+    // Recursive call
+    douglasPeuckerRecurse(list, epsilon, keep, begin, index);
+    douglasPeuckerRecurse(list, epsilon, keep, index, end);
+  }
+  else
+  {
+    // remove all points between begin and end
+    for (std::size_t i = begin + 1; i < end; ++i)
+      keep[i] = false;
+  }
+}
+
+FunctionNodeList
+douglasPeucker(const FunctionNodeList & list, const Real epsilon)
+{
+  // set up keep list for function nodes
+  std::vector<bool> keep(list.size(), true);
+  douglasPeuckerRecurse(list, epsilon, keep, 0, list.size() - 1);
+
+  FunctionNodeList result;
+  result.reserve(std::count_if(keep.begin(), keep.end(), [](bool k) { return k; }));
+
+  /// filter result
+  for (std::size_t i = 0; i < list.size(); ++i)
+    if (keep[i])
+      result.push_back(list[i]);
+
+  return result;
+}
+
+} // namespace PointReduction

framework/src/vectorpostprocessors/PiecewiseFunctionTabulate.C

@@ -0,0 +1,53 @@
+//* 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 "PiecewiseFunctionTabulate.h"
+#include "PiecewiseBase.h"
+
+registerMooseObject("MooseApp", PiecewiseFunctionTabulate);
+
+template <>
+InputParameters
+validParams<PiecewiseFunctionTabulate>()
+{
+  InputParameters params = validParams<GeneralVectorPostprocessor>();
+  params.addClassDescription("Tabulate the function nodes of a piecewise function, such as "
+                             "PiecewiseLinear or PiecewiseConstant");
+  params.addRequiredParam<FunctionName>(
+      "function", "Name of the piecewise function object to extract the time steps from");
+  return params;
+}
+
+PiecewiseFunctionTabulate::PiecewiseFunctionTabulate(const InputParameters & parameters)
+  : GeneralVectorPostprocessor(parameters),
+    _piecewise_function(dynamic_cast<PiecewiseBase *>(
+        &_fe_problem.getFunction(getParam<FunctionName>("function"),
+                                 isParamValid("_tid") ? getParam<THREAD_ID>("_tid") : 0))),
+    _x_col(declareVector("x")),
+    _y_col(declareVector("y"))
+{
+  if (!_piecewise_function)
+    paramError("function",
+               "The supplied function must be derived from PiecewiseBase (e.g. PiecewiseLinear or "
+               "PiecewiseConstant)");
+}
+
+void
+PiecewiseFunctionTabulate::execute()
+{
+  auto size = _piecewise_function->functionSize();
+  _x_col.resize(size);
+  _y_col.resize(size);
+
+  for (std::size_t i = 0; i < size; ++i)
+  {
+    _x_col[i] = _piecewise_function->domain(i);
+    _y_col[i] = _piecewise_function->range(i);
+  }
+}