From 454a08167dbc5d67dc43a6d91e1993e7b8405c07 Mon Sep 17 00:00:00 2001
From: Luca Heltai <luca.heltai@sissa.it>
Date: Thu, 18 Apr 2019 19:57:37 +0200
Subject: [PATCH] Added FunctionCutOffTensorProduct.

---
 include/deal.II/base/function_lib.h  | 109 +++++++++++++++++++++++-
 source/base/function_lib_cutoff.cc   |  79 ++++++++++++++++++
 tests/base/function_cutoff_02.cc     | 120 +++++++++++++++++++++++++++
 tests/base/function_cutoff_02.output |  64 ++++++++++++++
 4 files changed, 369 insertions(+), 3 deletions(-)
 create mode 100644 tests/base/function_cutoff_02.cc
 create mode 100644 tests/base/function_cutoff_02.output

diff --git a/include/deal.II/base/function_lib.h b/include/deal.II/base/function_lib.h
index 0b36b2d4e3..ff38a07bc8 100644
--- a/include/deal.II/base/function_lib.h
+++ b/include/deal.II/base/function_lib.h
@@ -21,6 +21,7 @@
 
 #include <deal.II/base/function.h>
 #include <deal.II/base/point.h>
+#include <deal.II/base/smartpointer.h>
 #include <deal.II/base/table.h>
 
 #include <array>
@@ -946,6 +947,11 @@ namespace Functions
       const bool         integrate_to_one       = false,
       const double       unitary_integral_value = 1.0);
 
+    /**
+     * Virtual destructor.
+     */
+    virtual ~CutOffFunctionBase() = default;
+
     /**
      * Move the center of the ball to new point <tt>p</tt>.
      *
@@ -965,7 +971,7 @@ namespace Functions
     /**
      * Set the center of the ball to the point @p p.
      */
-    void
+    virtual void
     set_center(const Point<dim> &p);
 
     /**
@@ -973,7 +979,7 @@ namespace Functions
      *
      * @deprecated Use set_radius() instead.
      */
-    void
+    virtual void
     set_radius(const double r);
 
     /**
@@ -988,6 +994,12 @@ namespace Functions
     double
     get_radius() const;
 
+    /**
+     * Return a boolean indicating if this function integrates to one.
+     */
+    bool
+    integrates_to_one() const;
+
   protected:
     /**
      * Center of the integration ball.
@@ -1008,7 +1020,7 @@ namespace Functions
     /**
      * Flag that controls wether we rescale the value when the radius changes.
      */
-    const bool integrate_to_one;
+    bool integrate_to_one;
 
     /**
      * The reference integral value. Derived classes should specify what their
@@ -1023,6 +1035,76 @@ namespace Functions
   };
 
 
+  /**
+   * Tensor product of CutOffFunctionBase objects.
+   *
+   * Instead of using the distance to compute the cut-off function, this class
+   * performs a tensor product of the same CutOffFunctionBase object in each
+   * coordinate direction.
+   *
+   * @ingroup functions
+   * @author Luca Heltai, 2019.
+   */
+  template <int dim>
+  class CutOffFunctionTensorProduct : public CutOffFunctionBase<dim>
+  {
+  public:
+    /**
+     * Construct an empty CutOffFunctionTensorProduct object.
+     *
+     * Before you can use this class, you have to call the set_base() method
+     * with a class derived from the CutOffFunctionBase object.
+     *
+     * If you try to use this class before you call the set_base() method,
+     * and exception will be triggered.
+     */
+    CutOffFunctionTensorProduct(
+      double             radius       = 1.0,
+      const Point<dim> & center       = Point<dim>(),
+      const unsigned int n_components = 1,
+      const unsigned int select       = CutOffFunctionBase<dim>::no_component,
+      const bool         integrate_to_one = false);
+
+    /**
+     * Initialize the class with an objet of type
+     * @tparam CutOffFunctionBaseType<1>.
+     */
+    template <template <int> class CutOffFunctionBaseType>
+    void
+    set_base();
+
+    /**
+     * Set the new center.
+     */
+    virtual void
+    set_center(const Point<dim> &center) override;
+
+    /**
+     * Set the new radius.
+     */
+    virtual void
+    set_radius(const double radius) override;
+
+    /**
+     * Function value at one point.
+     */
+    virtual double
+    value(const Point<dim> &p, const unsigned int component = 0) const override;
+
+    /**
+     * Function gradient at one point.
+     */
+    virtual Tensor<1, dim>
+    gradient(const Point<dim> & p,
+             const unsigned int component = 0) const override;
+
+  private:
+    std::array<std::unique_ptr<CutOffFunctionBase<1>>, dim> base;
+
+    bool initialized;
+  };
+
+
 
   /**
    * Cut-off function in L-infinity for an arbitrary ball.  This function is
@@ -1500,8 +1582,29 @@ namespace Functions
     const std::vector<double> coefficients;
   };
 
+#ifndef DOXYGEN
 
 
+
+  // Template definitions
+  template <int dim>
+  template <template <int> class CutOffFunctionBaseType>
+  void
+  CutOffFunctionTensorProduct<dim>::set_base()
+  {
+    initialized = true;
+    for (unsigned int i = 0; i < dim; ++i)
+      base[i].reset(new CutOffFunctionBaseType<1>(this->radius,
+                                                  Point<1>(this->center[i]),
+                                                  this->n_components,
+                                                  this->selected,
+                                                  this->integrate_to_one));
+  }
+
+
+
+#endif
+
 } // namespace Functions
 DEAL_II_NAMESPACE_CLOSE
 
diff --git a/source/base/function_lib_cutoff.cc b/source/base/function_lib_cutoff.cc
index e1fa60fc53..ff6ca15b0e 100644
--- a/source/base/function_lib_cutoff.cc
+++ b/source/base/function_lib_cutoff.cc
@@ -107,6 +107,85 @@ namespace Functions
   }
 
 
+
+  template <int dim>
+  bool
+  CutOffFunctionBase<dim>::integrates_to_one() const
+  {
+    return integrate_to_one;
+  }
+
+
+
+  template <int dim>
+  CutOffFunctionTensorProduct<dim>::CutOffFunctionTensorProduct(
+    double             radius,
+    const Point<dim> & center,
+    const unsigned int n_components,
+    const unsigned int select,
+    const bool         integrate_to_one)
+    : CutOffFunctionBase<dim>(radius,
+                              center,
+                              n_components,
+                              select,
+                              integrate_to_one)
+    , initialized(false)
+  {}
+
+
+
+  template <int dim>
+  void
+  CutOffFunctionTensorProduct<dim>::set_center(const Point<dim> &p)
+  {
+    Assert(initialized, ExcNotInitialized());
+    for (unsigned int i = 0; i < dim; ++i)
+      base[i]->set_center(Point<1>(p[i]));
+    CutOffFunctionBase<dim>::set_center(p);
+  }
+
+
+
+  template <int dim>
+  void
+  CutOffFunctionTensorProduct<dim>::set_radius(const double r)
+  {
+    Assert(initialized, ExcNotInitialized());
+    for (unsigned int i = 0; i < dim; ++i)
+      base[i]->set_radius(r);
+    CutOffFunctionBase<dim>::set_radius(r);
+  }
+
+
+
+  template <int dim>
+  double
+  CutOffFunctionTensorProduct<dim>::value(const Point<dim> & p,
+                                          const unsigned int component) const
+  {
+    Assert(initialized, ExcNotInitialized());
+    double ret = 1.0;
+    for (unsigned int i = 0; i < dim; ++i)
+      ret *= base[i]->value(Point<1>(p[i]), component);
+    return ret;
+  }
+
+
+
+  template <int dim>
+  Tensor<1, dim>
+  CutOffFunctionTensorProduct<dim>::gradient(const Point<dim> & p,
+                                             const unsigned int component) const
+  {
+    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];
+    return ret;
+  }
+
+
+
   //////////////////////////////////////////////////////////////////////
   namespace
   {
diff --git a/tests/base/function_cutoff_02.cc b/tests/base/function_cutoff_02.cc
new file mode 100644
index 0000000000..5540e56745
--- /dev/null
+++ b/tests/base/function_cutoff_02.cc
@@ -0,0 +1,120 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2019 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.
+//
+// ---------------------------------------------------------------------
+
+// Test the functionality of the CutOffTensorProductFunction
+
+#include <deal.II/base/function_lib.h>
+#include <deal.II/base/point.h>
+#include <deal.II/base/utilities.h>
+
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/tria.h>
+
+#include <deal.II/lac/vector.h>
+
+#include <deal.II/numerics/vector_tools.h>
+
+#include <map>
+
+#include "../tests.h"
+
+using namespace Functions;
+
+template <int dim>
+struct Domain
+{
+  Domain()
+    : fe(1)
+    , dh(tria)
+    , quad(5)
+  {
+    deallog << "Testing dim = " << dim << std::endl;
+    GridGenerator::hyper_cube(tria, -2, 2);
+    tria.refine_global(5 - dim);
+    dh.distribute_dofs(fe);
+    zero.reinit(dh.n_dofs());
+    cell_integral.reinit(tria.n_active_cells());
+  }
+
+  template <class Function>
+  void
+  integrate(const Function &fun)
+  {
+    deallog << "Integrating " << Utilities::type_to_string(fun)
+            << " -2,2 cube: ";
+
+    VectorTools::integrate_difference(
+      dh, zero, fun, cell_integral, quad, VectorTools::L1_norm);
+    const double integral =
+      VectorTools::compute_global_error(tria,
+                                        cell_integral,
+                                        VectorTools::L1_norm);
+
+    deallog << integral << std::endl;
+  }
+
+  Triangulation<dim> tria;
+  FE_Q<dim>          fe;
+  DoFHandler<dim>    dh;
+  Vector<double>     zero;
+  Vector<double>     cell_integral;
+  QGauss<dim>        quad;
+};
+
+template <int dim, template <int> class TestCutOffFunction>
+void
+test()
+{
+  static Domain<dim> domain;
+
+  CutOffFunctionTensorProduct<dim> fun(
+    1, Point<dim>(), 1, CutOffFunctionBase<dim>::no_component, true);
+  fun.template set_base<TestCutOffFunction>();
+
+  deallog << "Center: " << fun.get_center() << std::endl
+          << "Radius: " << fun.get_radius() << std::endl;
+
+  domain.integrate(fun);
+
+  Point<dim> new_center;
+  for (unsigned int i = 0; i < dim; ++i)
+    new_center[i] = .5;
+
+  fun.set_center(new_center);
+  fun.set_radius(.5);
+
+  deallog << "Center: " << fun.get_center() << std::endl
+          << "Radius: " << fun.get_radius() << std::endl;
+
+  domain.integrate(fun);
+}
+
+int
+main()
+{
+  initlog(1);
+
+  test<1, Functions::CutOffFunctionLinfty>();
+  test<2, Functions::CutOffFunctionLinfty>();
+  test<3, Functions::CutOffFunctionLinfty>();
+
+  test<1, Functions::CutOffFunctionW1>();
+  test<2, Functions::CutOffFunctionW1>();
+  test<3, Functions::CutOffFunctionW1>();
+
+  test<1, Functions::CutOffFunctionCinfty>();
+  test<2, Functions::CutOffFunctionCinfty>();
+  test<3, Functions::CutOffFunctionCinfty>();
+}
diff --git a/tests/base/function_cutoff_02.output b/tests/base/function_cutoff_02.output
new file mode 100644
index 0000000000..2f3fe3da81
--- /dev/null
+++ b/tests/base/function_cutoff_02.output
@@ -0,0 +1,64 @@
+
+DEAL::Testing dim = 1
+DEAL::Center: 0.00000
+DEAL::Radius: 1.00000
+DEAL::Integrating dealii::Functions::CutOffFunctionTensorProduct<1> -2,2 cube: 1.00000
+DEAL::Center: 0.500000
+DEAL::Radius: 0.500000
+DEAL::Integrating dealii::Functions::CutOffFunctionTensorProduct<1> -2,2 cube: 1.00000
+DEAL::Testing dim = 2
+DEAL::Center: 0.00000 0.00000
+DEAL::Radius: 1.00000
+DEAL::Integrating dealii::Functions::CutOffFunctionTensorProduct<2> -2,2 cube: 1.00000
+DEAL::Center: 0.500000 0.500000
+DEAL::Radius: 0.500000
+DEAL::Integrating dealii::Functions::CutOffFunctionTensorProduct<2> -2,2 cube: 1.00000
+DEAL::Testing dim = 3
+DEAL::Center: 0.00000 0.00000 0.00000
+DEAL::Radius: 1.00000
+DEAL::Integrating dealii::Functions::CutOffFunctionTensorProduct<3> -2,2 cube: 1.00000
+DEAL::Center: 0.500000 0.500000 0.500000
+DEAL::Radius: 0.500000
+DEAL::Integrating dealii::Functions::CutOffFunctionTensorProduct<3> -2,2 cube: 1.00000
+DEAL::Testing dim = 1
+DEAL::Center: 0.00000
+DEAL::Radius: 1.00000
+DEAL::Integrating dealii::Functions::CutOffFunctionTensorProduct<1> -2,2 cube: 1.00000
+DEAL::Center: 0.500000
+DEAL::Radius: 0.500000
+DEAL::Integrating dealii::Functions::CutOffFunctionTensorProduct<1> -2,2 cube: 1.00000
+DEAL::Testing dim = 2
+DEAL::Center: 0.00000 0.00000
+DEAL::Radius: 1.00000
+DEAL::Integrating dealii::Functions::CutOffFunctionTensorProduct<2> -2,2 cube: 1.00000
+DEAL::Center: 0.500000 0.500000
+DEAL::Radius: 0.500000
+DEAL::Integrating dealii::Functions::CutOffFunctionTensorProduct<2> -2,2 cube: 1.00000
+DEAL::Testing dim = 3
+DEAL::Center: 0.00000 0.00000 0.00000
+DEAL::Radius: 1.00000
+DEAL::Integrating dealii::Functions::CutOffFunctionTensorProduct<3> -2,2 cube: 1.00000
+DEAL::Center: 0.500000 0.500000 0.500000
+DEAL::Radius: 0.500000
+DEAL::Integrating dealii::Functions::CutOffFunctionTensorProduct<3> -2,2 cube: 1.17474
+DEAL::Testing dim = 1
+DEAL::Center: 0.00000
+DEAL::Radius: 1.00000
+DEAL::Integrating dealii::Functions::CutOffFunctionTensorProduct<1> -2,2 cube: 1.00002
+DEAL::Center: 0.500000
+DEAL::Radius: 0.500000
+DEAL::Integrating dealii::Functions::CutOffFunctionTensorProduct<1> -2,2 cube: 0.999735
+DEAL::Testing dim = 2
+DEAL::Center: 0.00000 0.00000
+DEAL::Radius: 1.00000
+DEAL::Integrating dealii::Functions::CutOffFunctionTensorProduct<2> -2,2 cube: 0.999471
+DEAL::Center: 0.500000 0.500000
+DEAL::Radius: 0.500000
+DEAL::Integrating dealii::Functions::CutOffFunctionTensorProduct<2> -2,2 cube: 1.00322
+DEAL::Testing dim = 3
+DEAL::Center: 0.00000 0.00000 0.00000
+DEAL::Radius: 1.00000
+DEAL::Integrating dealii::Functions::CutOffFunctionTensorProduct<3> -2,2 cube: 1.00484
+DEAL::Center: 0.500000 0.500000 0.500000
+DEAL::Radius: 0.500000
+DEAL::Integrating dealii::Functions::CutOffFunctionTensorProduct<3> -2,2 cube: 1.00779
-- 
2.39.5