From 772651e912cf3caf2394fddbef1cabb58ad2f832 Mon Sep 17 00:00:00 2001 From: Luca Heltai Date: Thu, 18 Apr 2019 09:31:50 +0200 Subject: [PATCH] Modify FunctionCutOff to support integration to one. --- doc/news/changes/minor/20190418LucaHeltai | 3 + include/deal.II/base/function_lib.h | 50 ++++++++- source/base/function_lib_cutoff.cc | 131 ++++++++++++++++++---- tests/function_cutoff_01.cc | 120 ++++++++++++++++++++ tests/function_cutoff_01.output | 64 +++++++++++ 5 files changed, 339 insertions(+), 29 deletions(-) create mode 100644 doc/news/changes/minor/20190418LucaHeltai create mode 100644 tests/function_cutoff_01.cc create mode 100644 tests/function_cutoff_01.output diff --git a/doc/news/changes/minor/20190418LucaHeltai b/doc/news/changes/minor/20190418LucaHeltai new file mode 100644 index 0000000000..a0b5e6a83b --- /dev/null +++ b/doc/news/changes/minor/20190418LucaHeltai @@ -0,0 +1,3 @@ +New: Functions::CutOffFunctionBase now supports rescaling the function to keep its integral equal to one. +
+(Luca Heltai, 2019/04/18) diff --git a/include/deal.II/base/function_lib.h b/include/deal.II/base/function_lib.h index 665cef4b3c..3769444c3e 100644 --- a/include/deal.II/base/function_lib.h +++ b/include/deal.II/base/function_lib.h @@ -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 class CutOffFunctionBase : public Function @@ -930,12 +934,17 @@ namespace Functions * * If an argument select 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 = Point(), const unsigned int n_components = 1, - const unsigned int select = CutOffFunctionBase::no_component); + const unsigned int select = CutOffFunctionBase::no_component, + const bool integrate_to_one = false, + const double unitary_integral_value = 1.0); /** * Move the center of the ball to new point p. @@ -949,6 +958,18 @@ namespace Functions void new_radius(const double r); + /** + * Return the center stored in this object. + */ + const Point & + 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 = Point(), const unsigned int n_components = 1, - const unsigned int select = CutOffFunctionBase::no_component); + const unsigned int select = CutOffFunctionBase::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 = Point(), const unsigned int n_components = 1, - const unsigned int select = CutOffFunctionBase::no_component); + const unsigned int select = CutOffFunctionBase::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 = Point(), const unsigned int n_components = 1, - const unsigned int select = CutOffFunctionBase::no_component); + const unsigned int select = CutOffFunctionBase::no_component, + bool integrate_to_one = false); /** * Function value at one point. diff --git a/source/base/function_lib_cutoff.cc b/source/base/function_lib_cutoff.cc index b35fcc42d9..71ba7b7834 100644 --- a/source/base/function_lib_cutoff.cc +++ b/source/base/function_lib_cutoff.cc @@ -27,15 +27,25 @@ DEAL_II_NAMESPACE_OPEN namespace Functions { template - CutOffFunctionBase::CutOffFunctionBase(const double r, - const Point p, - const unsigned int n_components, - const unsigned int select) + CutOffFunctionBase::CutOffFunctionBase( + const double r, + const Point p, + const unsigned int n_components, + const unsigned int select, + const bool integrate_to_one, + const double unitary_integral_value) : Function(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(radius)) : + 1.0) + { + Assert(r > 0, ExcMessage("You must specify a radius > 0.")); + } template @@ -46,22 +56,68 @@ namespace Functions } + + template + const Point & + CutOffFunctionBase::get_center() const + { + return center; + } + + + template void CutOffFunctionBase::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(radius)); + else + rescaling = 1.0; + } + + + + template + double + CutOffFunctionBase::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 CutOffFunctionLinfty::CutOffFunctionLinfty( const double r, const Point p, const unsigned int n_components, - const unsigned int select) - : CutOffFunctionBase(r, p, n_components, select) + const unsigned int select, + const bool integrate_to_one) + : CutOffFunctionBase(r, + p, + n_components, + select, + integrate_to_one, + integral_Linfty[dim - 1]) {} @@ -72,7 +128,7 @@ namespace Functions { if (this->selected == CutOffFunctionBase::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::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::no_component) values[k] = val; else @@ -121,13 +180,18 @@ namespace Functions } } - template CutOffFunctionW1::CutOffFunctionW1(const double r, const Point p, const unsigned int n_components, - const unsigned int select) - : CutOffFunctionBase(r, p, n_components, select) + const unsigned int select, + const bool integrate_to_one) + : CutOffFunctionBase(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::no_component) values[k] = val; else @@ -197,8 +268,14 @@ namespace Functions const double r, const Point p, const unsigned int n_components, - const unsigned int select) - : CutOffFunctionBase(r, p, n_components, select) + const unsigned int select, + bool integrate_to_one) + : CutOffFunctionBase(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::no_component) @@ -298,11 +376,16 @@ namespace Functions return ((e < -50) ? Point() : (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 index 0000000000..4286608630 --- /dev/null +++ b/tests/function_cutoff_01.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 CutOffFunction with integrate_to_one = true + +#include +#include +#include + +#include +#include + +#include + +#include + +#include + +#include "../tests.h" + +using namespace Functions; + +template +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 + 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 tria; + FE_Q fe; + DoFHandler dh; + Vector zero; + Vector cell_integral; + QGauss quad; +}; + +template class TestCutOffFunction> +void +test() +{ + static Domain domain; + TestCutOffFunction fun(1., + Point(), + 1, + Functions::CutOffFunctionBase::no_component, + /*integrate_to_one = */ true); + deallog << "Center: " << fun.get_center() << std::endl + << "Radius: " << fun.get_radius() << std::endl; + + domain.integrate(fun); + + Point 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 index 0000000000..2554c237b7 --- /dev/null +++ b/tests/function_cutoff_01.output @@ -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 -- 2.39.5