]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Modify FunctionCutOff to support integration to one.
authorLuca Heltai <luca.heltai@sissa.it>
Thu, 18 Apr 2019 07:31:50 +0000 (09:31 +0200)
committerLuca Heltai <luca.heltai@sissa.it>
Tue, 23 Apr 2019 12:59:00 +0000 (14:59 +0200)
doc/news/changes/minor/20190418LucaHeltai [new file with mode: 0644]
include/deal.II/base/function_lib.h
source/base/function_lib_cutoff.cc
tests/function_cutoff_01.cc [new file with mode: 0644]
tests/function_cutoff_01.output [new file with mode: 0644]

diff --git a/doc/news/changes/minor/20190418LucaHeltai b/doc/news/changes/minor/20190418LucaHeltai
new file mode 100644 (file)
index 0000000..a0b5e6a
--- /dev/null
@@ -0,0 +1,3 @@
+New: Functions::CutOffFunctionBase now supports rescaling the function to keep its integral equal to one. 
+<br>
+(Luca Heltai, 2019/04/18)
index 665cef4b3c8db4e4d3987e54aba66b77df78449e..3769444c3e746a3f986fed9c558129e5cc7dd153 100644 (file)
@@ -912,8 +912,12 @@ namespace Functions
    * radius of the supporting ball of a cut-off function. It also stores the
    * number of the non-zero component, if the function is vector-valued.
    *
+   * This class can also be used for approximated Dirac delta functions. These
+   * are special cut-off functions whose integral is always equal to one,
+   * independently on the radius of the supporting ball.
+   *
    * @ingroup functions
-   * @author Guido Kanschat, 2002
+   * @author Guido Kanschat, 2002, Luca Heltai, 2019.
    */
   template <int dim>
   class CutOffFunctionBase : public Function<dim>
@@ -930,12 +934,17 @@ namespace Functions
      *
      * If an argument <tt>select</tt> is given and not -1, the cut-off
      * function will be non-zero for this component only.
+     *
+     * If the argument @p integrate_to_one is set to true, then the value of
+     * the function is rescaled whenever a new radius is set.
      */
     CutOffFunctionBase(
       const double radius             = 1.,
       const Point<dim>                = Point<dim>(),
       const unsigned int n_components = 1,
-      const unsigned int select       = CutOffFunctionBase<dim>::no_component);
+      const unsigned int select       = CutOffFunctionBase<dim>::no_component,
+      const bool         integrate_to_one       = false,
+      const double       unitary_integral_value = 1.0);
 
     /**
      * Move the center of the ball to new point <tt>p</tt>.
@@ -949,6 +958,18 @@ namespace Functions
     void
     new_radius(const double r);
 
+    /**
+     * Return the center stored in this object.
+     */
+    const Point<dim> &
+    get_center() const;
+
+    /**
+     * Return the radius stored in this object.
+     */
+    double
+    get_radius() const;
+
   protected:
     /**
      * Center of the integration ball.
@@ -965,6 +986,22 @@ namespace Functions
      * in all components.
      */
     const unsigned int selected;
+
+    /**
+     * Flag that controls wether we rescale the value when the radius changes.
+     */
+    const bool integrate_to_one;
+
+    /**
+     * The reference integral value. Derived classes should specify what their
+     * integral is when @p radius = 1.0.
+     */
+    const double unitary_integral_value;
+
+    /**
+     * Current rescaling to apply the the cut-off function.
+     */
+    double rescaling;
   };
 
 
@@ -992,7 +1029,8 @@ namespace Functions
       const double radius             = 1.,
       const Point<dim>                = Point<dim>(),
       const unsigned int n_components = 1,
-      const unsigned int select       = CutOffFunctionBase<dim>::no_component);
+      const unsigned int select       = CutOffFunctionBase<dim>::no_component,
+      const bool         integrate_to_one = false);
 
     /**
      * Function value at one point.
@@ -1040,7 +1078,8 @@ namespace Functions
       const double radius             = 1.,
       const Point<dim>                = Point<dim>(),
       const unsigned int n_components = 1,
-      const unsigned int select       = CutOffFunctionBase<dim>::no_component);
+      const unsigned int select       = CutOffFunctionBase<dim>::no_component,
+      const bool         integrate_to_one = false);
 
     /**
      * Function value at one point.
@@ -1089,7 +1128,8 @@ namespace Functions
       const double radius             = 1.,
       const Point<dim>                = Point<dim>(),
       const unsigned int n_components = 1,
-      const unsigned int select       = CutOffFunctionBase<dim>::no_component);
+      const unsigned int select       = CutOffFunctionBase<dim>::no_component,
+      bool               integrate_to_one = false);
 
     /**
      * Function value at one point.
index b35fcc42d9ecd6505552eebc4da72baae91f4274..71ba7b7834fc1f7e703c3232ae66939063270718 100644 (file)
@@ -27,15 +27,25 @@ DEAL_II_NAMESPACE_OPEN
 namespace Functions
 {
   template <int dim>
-  CutOffFunctionBase<dim>::CutOffFunctionBase(const double       r,
-                                              const Point<dim>   p,
-                                              const unsigned int n_components,
-                                              const unsigned int select)
+  CutOffFunctionBase<dim>::CutOffFunctionBase(
+    const double       r,
+    const Point<dim>   p,
+    const unsigned int n_components,
+    const unsigned int select,
+    const bool         integrate_to_one,
+    const double       unitary_integral_value)
     : Function<dim>(n_components)
     , center(p)
     , radius(r)
     , selected(select)
-  {}
+    , integrate_to_one(integrate_to_one)
+    , unitary_integral_value(unitary_integral_value)
+    , rescaling(integrate_to_one ? 1. / (unitary_integral_value *
+                                         Utilities::fixed_power<dim>(radius)) :
+                                   1.0)
+  {
+    Assert(r > 0, ExcMessage("You must specify a radius > 0."));
+  }
 
 
   template <int dim>
@@ -46,22 +56,68 @@ namespace Functions
   }
 
 
+
+  template <int dim>
+  const Point<dim> &
+  CutOffFunctionBase<dim>::get_center() const
+  {
+    return center;
+  }
+
+
+
   template <int dim>
   void
   CutOffFunctionBase<dim>::new_radius(const double r)
   {
     radius = r;
+    Assert(r > 0, ExcMessage("You must specify a radius > 0."));
+    if (integrate_to_one)
+      rescaling =
+        1. / (unitary_integral_value * Utilities::fixed_power<dim>(radius));
+    else
+      rescaling = 1.0;
+  }
+
+
+
+  template <int dim>
+  double
+  CutOffFunctionBase<dim>::get_radius() const
+  {
+    return radius;
   }
 
+
   //////////////////////////////////////////////////////////////////////
+  namespace
+  {
+    double integral_Linfty[] = {2.0,
+                                3.14159265358979323846264338328,
+                                4.18879020478639098461685784437};
+    double integral_W1[]     = {1.0,
+                            1.04719755119659774615421446109,
+                            1.04719755119659774615421446109};
+
+    double integral_Cinfty[] = {1.20690032243787617533623799633,
+                                1.26811216112759608094632335664,
+                                1.1990039070192139033798473858};
+  } // namespace
+
 
   template <int dim>
   CutOffFunctionLinfty<dim>::CutOffFunctionLinfty(
     const double       r,
     const Point<dim>   p,
     const unsigned int n_components,
-    const unsigned int select)
-    : CutOffFunctionBase<dim>(r, p, n_components, select)
+    const unsigned int select,
+    const bool         integrate_to_one)
+    : CutOffFunctionBase<dim>(r,
+                              p,
+                              n_components,
+                              select,
+                              integrate_to_one,
+                              integral_Linfty[dim - 1])
   {}
 
 
@@ -72,7 +128,7 @@ namespace Functions
   {
     if (this->selected == CutOffFunctionBase<dim>::no_component ||
         component == this->selected)
-      return ((this->center.distance(p) < this->radius) ? 1. : 0.);
+      return ((this->center.distance(p) < this->radius) ? this->rescaling : 0.);
     return 0.;
   }
 
@@ -92,7 +148,9 @@ namespace Functions
     if (this->selected == CutOffFunctionBase<dim>::no_component ||
         component == this->selected)
       for (unsigned int k = 0; k < values.size(); ++k)
-        values[k] = (this->center.distance(points[k]) < this->radius) ? 1. : 0.;
+        values[k] = (this->center.distance(points[k]) < this->radius) ?
+                      this->rescaling :
+                      0.;
     else
       std::fill(values.begin(), values.end(), 0.);
   }
@@ -109,8 +167,9 @@ namespace Functions
 
     for (unsigned int k = 0; k < values.size(); ++k)
       {
-        const double val =
-          (this->center.distance(points[k]) < this->radius) ? 1. : 0.;
+        const double val = (this->center.distance(points[k]) < this->radius) ?
+                             this->rescaling :
+                             0.;
         if (this->selected == CutOffFunctionBase<dim>::no_component)
           values[k] = val;
         else
@@ -121,13 +180,18 @@ namespace Functions
       }
   }
 
-
   template <int dim>
   CutOffFunctionW1<dim>::CutOffFunctionW1(const double       r,
                                           const Point<dim>   p,
                                           const unsigned int n_components,
-                                          const unsigned int select)
-    : CutOffFunctionBase<dim>(r, p, n_components, select)
+                                          const unsigned int select,
+                                          const bool         integrate_to_one)
+    : CutOffFunctionBase<dim>(r,
+                              p,
+                              n_components,
+                              select,
+                              integrate_to_one,
+                              integral_W1[dim - 1])
   {}
 
 
@@ -140,7 +204,9 @@ namespace Functions
         component == this->selected)
       {
         const double d = this->center.distance(p);
-        return ((d < this->radius) ? (this->radius - d) : 0.);
+        return ((d < this->radius) ?
+                  (this->radius - d) / this->radius * this->rescaling :
+                  0.);
       }
     return 0.;
   }
@@ -160,7 +226,9 @@ namespace Functions
       for (unsigned int i = 0; i < values.size(); ++i)
         {
           const double d = this->center.distance(points[i]);
-          values[i]      = ((d < this->radius) ? (this->radius - d) : 0.);
+          values[i]      = ((d < this->radius) ?
+                         (this->radius - d) / this->radius * this->rescaling :
+                         0.);
         }
     else
       std::fill(values.begin(), values.end(), 0.);
@@ -179,8 +247,11 @@ namespace Functions
 
     for (unsigned int k = 0; k < values.size(); ++k)
       {
-        const double d   = this->center.distance(points[k]);
-        const double val = (d < this->radius) ? (this->radius - d) : 0.;
+        const double d = this->center.distance(points[k]);
+        const double val =
+          (d < this->radius) ?
+            (this->radius - d) / this->radius * this->rescaling :
+            0.;
         if (this->selected == CutOffFunctionBase<dim>::no_component)
           values[k] = val;
         else
@@ -197,8 +268,14 @@ namespace Functions
     const double       r,
     const Point<dim>   p,
     const unsigned int n_components,
-    const unsigned int select)
-    : CutOffFunctionBase<dim>(r, p, n_components, select)
+    const unsigned int select,
+    bool               integrate_to_one)
+    : CutOffFunctionBase<dim>(r,
+                              p,
+                              n_components,
+                              select,
+                              integrate_to_one,
+                              integral_Cinfty[dim - 1])
   {}
 
 
@@ -215,7 +292,7 @@ namespace Functions
         if (d >= r)
           return 0.;
         const double e = -r * r / (r * r - d * d);
-        return ((e < -50) ? 0. : numbers::E * std::exp(e));
+        return ((e < -50) ? 0. : numbers::E * std::exp(e) * this->rescaling);
       }
     return 0.;
   }
@@ -244,7 +321,8 @@ namespace Functions
           else
             {
               const double e = -r * r / (r * r - d * d);
-              values[i]      = (e < -50) ? 0. : numbers::E * std::exp(e);
+              values[i] =
+                (e < -50) ? 0. : numbers::E * std::exp(e) * this->rescaling;
             }
         }
     else
@@ -270,7 +348,7 @@ namespace Functions
           {
             const double e = -r * r / (r * r - d * d);
             if (e > -50)
-              val = numbers::E * std::exp(e);
+              val = numbers::E * std::exp(e) * this->rescaling;
           }
 
         if (this->selected == CutOffFunctionBase<dim>::no_component)
@@ -298,11 +376,16 @@ namespace Functions
     return ((e < -50) ? Point<dim>() :
                         (p - this->center) / d *
                           (-2.0 * r * r / std::pow(-r * r + d * d, 2.0) * d *
-                           std::exp(e)));
+                           std::exp(e)) *
+                          this->rescaling);
   }
 
 
   // explicit instantiations
+  template class CutOffFunctionBase<1>;
+  template class CutOffFunctionBase<2>;
+  template class CutOffFunctionBase<3>;
+
   template class CutOffFunctionLinfty<1>;
   template class CutOffFunctionLinfty<2>;
   template class CutOffFunctionLinfty<3>;
diff --git a/tests/function_cutoff_01.cc b/tests/function_cutoff_01.cc
new file mode 100644 (file)
index 0000000..4286608
--- /dev/null
@@ -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 CutOffFunction with integrate_to_one = true
+
+#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;
+  TestCutOffFunction<dim> fun(1.,
+                              Point<dim>(),
+                              1,
+                              Functions::CutOffFunctionBase<dim>::no_component,
+                              /*integrate_to_one = */ true);
+  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.new_center(new_center);
+  fun.new_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/function_cutoff_01.output b/tests/function_cutoff_01.output
new file mode 100644 (file)
index 0000000..2554c23
--- /dev/null
@@ -0,0 +1,64 @@
+
+DEAL::Testing dim = 1
+DEAL::Center: 0.00000
+DEAL::Radius: 1.00000
+DEAL::Integrating dealii::Functions::CutOffFunctionLinfty<1> -2,2 cube: 1.00000
+DEAL::Center: 0.500000
+DEAL::Radius: 0.500000
+DEAL::Integrating dealii::Functions::CutOffFunctionLinfty<1> -2,2 cube: 1.00000
+DEAL::Testing dim = 2
+DEAL::Center: 0.00000 0.00000
+DEAL::Radius: 1.00000
+DEAL::Integrating dealii::Functions::CutOffFunctionLinfty<2> -2,2 cube: 0.993980
+DEAL::Center: 0.500000 0.500000
+DEAL::Radius: 0.500000
+DEAL::Integrating dealii::Functions::CutOffFunctionLinfty<2> -2,2 cube: 1.02445
+DEAL::Testing dim = 3
+DEAL::Center: 0.00000 0.00000 0.00000
+DEAL::Radius: 1.00000
+DEAL::Integrating dealii::Functions::CutOffFunctionLinfty<3> -2,2 cube: 1.00850
+DEAL::Center: 0.500000 0.500000 0.500000
+DEAL::Radius: 0.500000
+DEAL::Integrating dealii::Functions::CutOffFunctionLinfty<3> -2,2 cube: 0.958427
+DEAL::Testing dim = 1
+DEAL::Center: 0.00000
+DEAL::Radius: 1.00000
+DEAL::Integrating dealii::Functions::CutOffFunctionW1<1> -2,2 cube: 1.00000
+DEAL::Center: 0.500000
+DEAL::Radius: 0.500000
+DEAL::Integrating dealii::Functions::CutOffFunctionW1<1> -2,2 cube: 1.00000
+DEAL::Testing dim = 2
+DEAL::Center: 0.00000 0.00000
+DEAL::Radius: 1.00000
+DEAL::Integrating dealii::Functions::CutOffFunctionW1<2> -2,2 cube: 1.00030
+DEAL::Center: 0.500000 0.500000
+DEAL::Radius: 0.500000
+DEAL::Integrating dealii::Functions::CutOffFunctionW1<2> -2,2 cube: 1.00316
+DEAL::Testing dim = 3
+DEAL::Center: 0.00000 0.00000 0.00000
+DEAL::Radius: 1.00000
+DEAL::Integrating dealii::Functions::CutOffFunctionW1<3> -2,2 cube: 1.00445
+DEAL::Center: 0.500000 0.500000 0.500000
+DEAL::Radius: 0.500000
+DEAL::Integrating dealii::Functions::CutOffFunctionW1<3> -2,2 cube: 1.03922
+DEAL::Testing dim = 1
+DEAL::Center: 0.00000
+DEAL::Radius: 1.00000
+DEAL::Integrating dealii::Functions::CutOffFunctionCinfty<1> -2,2 cube: 1.00002
+DEAL::Center: 0.500000
+DEAL::Radius: 0.500000
+DEAL::Integrating dealii::Functions::CutOffFunctionCinfty<1> -2,2 cube: 0.999735
+DEAL::Testing dim = 2
+DEAL::Center: 0.00000 0.00000
+DEAL::Radius: 1.00000
+DEAL::Integrating dealii::Functions::CutOffFunctionCinfty<2> -2,2 cube: 0.999958
+DEAL::Center: 0.500000 0.500000
+DEAL::Radius: 0.500000
+DEAL::Integrating dealii::Functions::CutOffFunctionCinfty<2> -2,2 cube: 0.998842
+DEAL::Testing dim = 3
+DEAL::Center: 0.00000 0.00000 0.00000
+DEAL::Radius: 1.00000
+DEAL::Integrating dealii::Functions::CutOffFunctionCinfty<3> -2,2 cube: 0.999614
+DEAL::Center: 0.500000 0.500000 0.500000
+DEAL::Radius: 0.500000
+DEAL::Integrating dealii::Functions::CutOffFunctionCinfty<3> -2,2 cube: 1.00165

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.