From e0056431e6be413a0fc91d095ac5f6651f126900 Mon Sep 17 00:00:00 2001
From: Sebastian Schunert <snschune@ncsu.edu>
Date: Tue, 7 Jan 2020 18:12:44 -0700
Subject: [PATCH] Bug fix and test (#410)

---
 .../PolyatomicDisplacementFunctionBase.C      |  4 +-
 unit/src/DisplacementFunctionTest.C           | 59 +++++++++++++++++++
 2 files changed, 62 insertions(+), 1 deletion(-)
 create mode 100644 unit/src/DisplacementFunctionTest.C

diff --git a/src/utils/PolyatomicDisplacementFunctionBase.C b/src/utils/PolyatomicDisplacementFunctionBase.C
index f6e6d86d..672a0934 100644
--- a/src/utils/PolyatomicDisplacementFunctionBase.C
+++ b/src/utils/PolyatomicDisplacementFunctionBase.C
@@ -126,7 +126,8 @@ PolyatomicDisplacementFunctionBase::computeDisplacementFunctionIntegral()
 
   for (unsigned int e = 1; e < nEnergySteps(); ++e)
   {
-    _displacement_function_integral.push_back(std::vector<Real>(_problem_size));
+    // initialize the integral to the integral of the previous step
+    _displacement_function_integral[e] = _displacement_function_integral[e - 1];
 
     // set energy points for integration
     Real lower = energyPoint(e - 1);
@@ -141,6 +142,7 @@ PolyatomicDisplacementFunctionBase::computeDisplacementFunctionIntegral()
       {
         unsigned int i, j, l;
         inverseMapIndex(n, i, j, l);
+
         _displacement_function_integral[e][n] += w * linearInterpolation(energy, i, j, l);
       }
     }
diff --git a/unit/src/DisplacementFunctionTest.C b/unit/src/DisplacementFunctionTest.C
new file mode 100644
index 00000000..416db1a4
--- /dev/null
+++ b/unit/src/DisplacementFunctionTest.C
@@ -0,0 +1,59 @@
+/****************************************************************/
+/*               DO NOT MODIFY THIS HEADER                      */
+/* MOOSE - Multiphysics Object Oriented Simulation Environment  */
+/*                                                              */
+/*           (c) 2010 Battelle Energy Alliance, LLC             */
+/*                   ALL RIGHTS RESERVED                        */
+/*                                                              */
+/*          Prepared by Battelle Energy Alliance, LLC           */
+/*            Under Contract No. DE-AC07-05ID14517              */
+/*            With the U. S. Department of Energy               */
+/*                                                              */
+/*            See COPYRIGHT for full restrictions               */
+/****************************************************************/
+
+#include "PolyatomicDisplacementFunction.h"
+#include <gtest/gtest.h>
+#include <mytrim/element.h>
+
+TEST(DisplacementFunctionTest, integralDispFunction)
+{
+  std::vector<unsigned int> Z = {6, 14};
+  std::vector<Real> A = {12.0, 28.0};
+  std::vector<Real> N = {0.5, 0.5};
+  std::vector<Real> threshold = {16.3, 92.6};
+  std::vector<MyTRIM_NS::Element> poly_mat;
+  for (unsigned int j = 0; j < Z.size(); ++j)
+  {
+    MyTRIM_NS::Element element;
+    element._Z = Z[j];
+    element._m = A[j];
+    element._t = N[j];
+    element._Edisp = threshold[j];
+    element._Elbind = 0.0;
+    poly_mat.push_back(element);
+  }
+
+  PolyatomicDisplacementFunction padf(poly_mat, NET);
+
+  Real energy = padf.minEnergy();
+  for (unsigned int j = 0; j < 51; ++j)
+  {
+    energy += 1;
+    padf.advanceDisplacements(energy);
+  }
+  padf.computeDisplacementFunctionIntegral();
+  unsigned int final = padf.nEnergySteps() - 1;
+
+  Real integral = 0;
+  for (unsigned int j = 0; j <= final; ++j)
+  {
+    Real mp = (j == 0 || j == final) ? 0.5 : 1;
+    Real energy = padf.energyPoint(j);
+    integral += mp * padf.linearInterpolation(energy, 0, 0, 0);
+  }
+  EXPECT_NEAR(integral,
+              padf.linearInterpolationIntegralDamageFunction(padf.energyPoint(final), 0, 0, 0),
+              1e-4)
+      << "Displacement function integration result is wrong";
+}