Fixed gradient of tensor product cutoff.

diff --git a/doc/news/changes/minor/20210213LucaHeltai b/doc/news/changes/minor/20210213LucaHeltai
new file mode 100644 (file)
index 0000000..bb1e472
--- /dev/null
+++ b/doc/news/changes/minor/20210213LucaHeltai
@@ -0,0 +1,3 @@
+Fixed: CutOffFunctionTensorProduct::gradient() was computing the wrong gradient. This is now fixed.
+<br>
+(Luca Heltai, 2021/02/13)

diff --git a/source/base/function_lib_cutoff.cc b/source/base/function_lib_cutoff.cc
index c3fbdb8507d12a3645b061638a65e86ae6ba01bf..03922d8d46c1f99188d43873124e4420bd1a57b5 100644 (file)
--- a/source/base/function_lib_cutoff.cc
+++ b/source/base/function_lib_cutoff.cc
@@ -162,8 +162,13 @@ namespace Functions
   {
     Assert(initialized, ExcNotInitialized());
     Tensor<1, dim> ret;
-    for (unsigned int i = 0; i < dim; ++i)
-      ret[i] = base[i]->gradient(Point<1>(p[i]), component)[0];
+    for (unsigned int d = 0; d < dim; ++d)
+      {
+        ret[d] = base[d]->gradient(Point<1>(p[d]), component)[0];
+        for (unsigned int i = 0; i < dim; ++i)
+          if (i != d)
+            ret[d] *= base[i]->value(Point<1>(p[i]), component);
+      }
     return ret;
   }
 

diff --git a/tests/base/function_cutoff_03.cc b/tests/base/function_cutoff_03.cc
new file mode 100644 (file)
index 0000000..7c1d3a9
--- /dev/null
+++ b/tests/base/function_cutoff_03.cc
@@ -0,0 +1,85 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2019 - 2020 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE.md at
+// the top level directory of deal.II.
+//
+// ---------------------------------------------------------------------
+
+// Check that CutOffTensorProductFunction produces the right gradient in dim>1
+
+#include <deal.II/base/function_lib.h>
+#include <deal.II/base/point.h>
+#include <deal.II/base/quadrature_lib.h>
+
+#include "../tests.h"
+
+using namespace Functions;
+
+template <int dim>
+void
+test()
+{
+  Point<dim> p;
+  for (unsigned int i = 0; i < dim; ++i)
+    p[i] = .5;
+
+
+  CutOffFunctionTensorProduct<dim> fun(
+    .5, p, 1, CutOffFunctionBase<dim>::no_component, true);
+
+
+  fun.template set_base<CutOffFunctionC1>();
+
+  deallog << "Center: " << fun.get_center() << std::endl
+          << "Radius: " << fun.get_radius() << std::endl;
+
+  // Check integration of the function, and of its gradient square
+  QGauss<dim> quad(12);
+
+  std::vector<Tensor<1, dim>> gradients(quad.size());
+  std::vector<double>         values(quad.size());
+  fun.value_list(quad.get_points(), values);
+  fun.gradient_list(quad.get_points(), gradients);
+
+  // Integral in 1d of grad square: 2.0*pi^2
+  // Integral in 2d of grad square: 6.0*pi^2
+  // Integral in 3d of grad square: 13.5*pi^2
+
+  const double spi = numbers::PI * numbers::PI;
+
+  const std::array<double, 3> exact{{2.0 * spi, 6.0 * spi, 13.5 * spi}};
+
+  double integral      = 0.0;
+  double grad_integral = 0.0;
+  for (unsigned int i = 0; i < quad.size(); ++i)
+    {
+      integral += values[i] * quad.weight(i);
+      grad_integral += gradients[i] * gradients[i] * quad.weight(i);
+    }
+  if (std::abs(integral - 1.0) > 1e-10)
+    deallog << "Value NOT OK" << std::endl;
+  else
+    deallog << "Value OK" << std::endl;
+  if (std::abs(grad_integral - exact[dim - 1]) > 1e-10)
+    deallog << "Gradient NOT OK" << std::endl;
+  else
+    deallog << "Gradient OK" << std::endl;
+}
+
+int
+main()
+{
+  initlog(1);
+
+  test<1>();
+  test<2>();
+  test<3>();
+}

diff --git a/tests/base/function_cutoff_03.output b/tests/base/function_cutoff_03.output
new file mode 100644 (file)
index 0000000..8c8ba15
--- /dev/null
+++ b/tests/base/function_cutoff_03.output
@@ -0,0 +1,13 @@
+
+DEAL::Center: 0.500000
+DEAL::Radius: 0.500000
+DEAL::Value OK
+DEAL::Gradient OK
+DEAL::Center: 0.500000 0.500000
+DEAL::Radius: 0.500000
+DEAL::Value OK
+DEAL::Gradient OK
+DEAL::Center: 0.500000 0.500000 0.500000
+DEAL::Radius: 0.500000
+DEAL::Value OK
+DEAL::Gradient OK