--- /dev/null
+New: Functions::CutOffFunctionBase now supports rescaling the function to keep its integral equal to one.
+<br>
+(Luca Heltai, 2019/04/18)
* 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>
*
* 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>.
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.
* 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;
};
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.
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.
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.
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>
}
+
+ 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])
{}
{
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.;
}
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.);
}
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
}
}
-
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])
{}
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.;
}
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.);
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
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])
{}
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.;
}
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
{
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)
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>;
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// 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>();
+}
--- /dev/null
+
+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
projects / dealii.git / commitdiff