From 142d7ab2210dd99c1e3056988ae52c06f747f02e Mon Sep 17 00:00:00 2001 From: Luca Heltai Date: Thu, 18 Apr 2019 21:42:08 +0200 Subject: [PATCH] Added implementation of CutOffFunctionC1 --- include/deal.II/base/function_lib.h | 85 ++++++++++++--- source/base/function_lib_cutoff.cc | 150 +++++++++++++++++++++++++-- tests/base/function_cutoff_01.cc | 4 + tests/base/function_cutoff_01.output | 21 ++++ tests/base/function_cutoff_02.cc | 4 + tests/base/function_cutoff_02.output | 21 ++++ 6 files changed, 263 insertions(+), 22 deletions(-) diff --git a/include/deal.II/base/function_lib.h b/include/deal.II/base/function_lib.h index ff38a07bc8..ab7635e89a 100644 --- a/include/deal.II/base/function_lib.h +++ b/include/deal.II/base/function_lib.h @@ -915,7 +915,7 @@ namespace Functions * * 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. + * independently of the radius of the supporting ball. * * @ingroup functions * @author Guido Kanschat, 2002, Luca Heltai, 2019. @@ -931,17 +931,23 @@ namespace Functions static const unsigned int no_component = numbers::invalid_unsigned_int; /** - * Constructor. Arguments are the center of the ball and its radius. - * - * If an argument select is given and not -1, the cut-off - * function will be non-zero for this component only. + * Constructor. * - * 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. + * @param[in] radius Radius of the ball + * @param[in] center Center of the ball + * @param[in] n_components Number of components of this function object + * @param[in] select If this is different from + * CutOffFunctionBase::no_component, then the function will be non-zero + * for this component only + * @param[in] integrate_to_one Rescale the value of the function whenever a + * new radius is set, to guarantee that the integral is equal to one + * @param[in] unitary_integral_value Value of the integral when the radius + * is equal to 1.0. Derived classes will need to supply this value, to + * guarantee that the rescaling is performed correctly. */ CutOffFunctionBase( - const double radius = 1., - const Point = Point(), + const double radius = 1., + const Point center = Point(), const unsigned int n_components = 1, const unsigned int select = CutOffFunctionBase::no_component, const bool integrate_to_one = false, @@ -976,8 +982,6 @@ namespace Functions /** * Set the radius of the ball to @p r - * - * @deprecated Use set_radius() instead. */ virtual void set_radius(const double r); @@ -995,7 +999,7 @@ namespace Functions get_radius() const; /** - * Return a boolean indicating if this function integrates to one. + * Return a boolean indicating whether this function integrates to one. */ bool integrates_to_one() const; @@ -1204,6 +1208,63 @@ namespace Functions }; + /** + * A cut-off function for an arbitrarily-sized ball that is in the space $C^1$ + * (i.e., continuously differentiable). This is a cut-off function that is + * often used in the literature of the Immersed Boundary Method. + * + * The expression of the function in radial coordinates is given by + * $f(r)=1/2(cos(\pi r/s)+1)$ where $r + class CutOffFunctionC1 : public CutOffFunctionBase + { + public: + /** + * Constructor. + */ + CutOffFunctionC1( + const double radius = 1., + const Point = Point(), + const unsigned int n_components = 1, + const unsigned int select = CutOffFunctionBase::no_component, + bool integrate_to_one = false); + + /** + * Function value at one point. + */ + virtual double + value(const Point &p, const unsigned int component = 0) const override; + + /** + * Function values at multiple points. + */ + virtual void + value_list(const std::vector> &points, + std::vector & values, + const unsigned int component = 0) const override; + + /** + * Function values at multiple points. + */ + virtual void + vector_value_list(const std::vector> &points, + std::vector> & values) const override; + + /** + * Function gradient at one point. + */ + virtual Tensor<1, dim> + gradient(const Point & p, + const unsigned int component = 0) const override; + }; + + /** * Cut-off function for an arbitrary ball. This is the traditional cut-off * function in C-infinity for a ball of certain radius around diff --git a/source/base/function_lib_cutoff.cc b/source/base/function_lib_cutoff.cc index ff6ca15b0e..87d79d9195 100644 --- a/source/base/function_lib_cutoff.cc +++ b/source/base/function_lib_cutoff.cc @@ -189,16 +189,29 @@ namespace Functions ////////////////////////////////////////////////////////////////////// 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}; + // Integral of CutOffFunctionLinfty in dimension 1, 2, and 3 when the radius + // is one + const double integral_Linfty[] = {2.0, + 3.14159265358979323846264338328, + 4.18879020478639098461685784437}; + + // Integral of CutOffFunctionW1 in dimension 1, 2, and 3 when the radius + // is one + const double integral_W1[] = {1.0, + 1.04719755119659774615421446109, + 1.04719755119659774615421446109}; + + // Integral of CutOffFunctionCinfty in dimension 1, 2, and 3 when the radius + // is one + const double integral_Cinfty[] = {1.20690032243787617533623799633, + 1.26811216112759608094632335664, + 1.1990039070192139033798473858}; + + // Integral of CutOffFunctionC1 in dimension 1, 2, and 3 when the radius + // is one + const double integral_C1[] = {1.0, + 0.93417655442731527615578663815, + 0.821155557658032806157358815206}; } // namespace @@ -478,6 +491,115 @@ namespace Functions } + + template + CutOffFunctionC1::CutOffFunctionC1(const double r, + const Point p, + const unsigned int n_components, + const unsigned int select, + bool integrate_to_one) + : CutOffFunctionBase(r, + p, + n_components, + select, + integrate_to_one, + integral_C1[dim - 1]) + {} + + + template + double + CutOffFunctionC1::value(const Point & p, + const unsigned int component) const + { + if (this->selected == CutOffFunctionBase::no_component || + component == this->selected) + { + const double d = this->center.distance(p); + const double r = this->radius; + if (d >= r) + return 0.; + return .5 * (std::cos(numbers::PI * d / r) + 1) * this->rescaling; + } + return 0.; + } + + + template + void + CutOffFunctionC1::value_list(const std::vector> &points, + std::vector & values, + const unsigned int component) const + { + Assert(values.size() == points.size(), + ExcDimensionMismatch(values.size(), points.size())); + + const double r = this->radius; + + if (this->selected == CutOffFunctionBase::no_component || + component == this->selected) + for (unsigned int i = 0; i < values.size(); ++i) + { + const double d = this->center.distance(points[i]); + if (d >= r) + { + values[i] = 0.; + } + else + { + values[i] = + .5 * (std::cos(numbers::PI * d / r) + 1) * this->rescaling; + } + } + else + std::fill(values.begin(), values.end(), 0.); + } + + + template + void + CutOffFunctionC1::vector_value_list( + const std::vector> &points, + std::vector> & values) const + { + Assert(values.size() == points.size(), + ExcDimensionMismatch(values.size(), points.size())); + + for (unsigned int k = 0; k < values.size(); ++k) + { + const double d = this->center.distance(points[k]); + const double r = this->radius; + double val = 0.; + if (d < this->radius) + { + val = .5 * (std::cos(numbers::PI * d / r) + 1) * this->rescaling; + } + + if (this->selected == CutOffFunctionBase::no_component) + values[k] = val; + else + { + values[k] = 0; + values[k](this->selected) = val; + } + } + } + + + + template + Tensor<1, dim> + CutOffFunctionC1::gradient(const Point &p, const unsigned int) const + { + const double d = this->center.distance(p); + const double r = this->radius; + if (d >= r) + return Tensor<1, dim>(); + return (-0.5 * numbers::PI * std::sin(numbers::PI * d / r) / r) * + (p - this->center) / d * this->rescaling; + } + + // explicit instantiations template class CutOffFunctionBase<1>; template class CutOffFunctionBase<2>; @@ -494,6 +616,14 @@ namespace Functions template class CutOffFunctionCinfty<1>; template class CutOffFunctionCinfty<2>; template class CutOffFunctionCinfty<3>; + + template class CutOffFunctionC1<1>; + template class CutOffFunctionC1<2>; + template class CutOffFunctionC1<3>; + + template class CutOffFunctionTensorProduct<1>; + template class CutOffFunctionTensorProduct<2>; + template class CutOffFunctionTensorProduct<3>; } // namespace Functions DEAL_II_NAMESPACE_CLOSE diff --git a/tests/base/function_cutoff_01.cc b/tests/base/function_cutoff_01.cc index c51ac0f797..2757d5bcf8 100644 --- a/tests/base/function_cutoff_01.cc +++ b/tests/base/function_cutoff_01.cc @@ -117,4 +117,8 @@ main() test<1, Functions::CutOffFunctionCinfty>(); test<2, Functions::CutOffFunctionCinfty>(); test<3, Functions::CutOffFunctionCinfty>(); + + test<1, Functions::CutOffFunctionC1>(); + test<2, Functions::CutOffFunctionC1>(); + test<3, Functions::CutOffFunctionC1>(); } diff --git a/tests/base/function_cutoff_01.output b/tests/base/function_cutoff_01.output index 2554c237b7..f36ce814d2 100644 --- a/tests/base/function_cutoff_01.output +++ b/tests/base/function_cutoff_01.output @@ -62,3 +62,24 @@ 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 +DEAL::Testing dim = 1 +DEAL::Center: 0.00000 +DEAL::Radius: 1.00000 +DEAL::Integrating dealii::Functions::CutOffFunctionC1<1> -2,2 cube: 1.00000 +DEAL::Center: 0.500000 +DEAL::Radius: 0.500000 +DEAL::Integrating dealii::Functions::CutOffFunctionC1<1> -2,2 cube: 1.00000 +DEAL::Testing dim = 2 +DEAL::Center: 0.00000 0.00000 +DEAL::Radius: 1.00000 +DEAL::Integrating dealii::Functions::CutOffFunctionC1<2> -2,2 cube: 0.999972 +DEAL::Center: 0.500000 0.500000 +DEAL::Radius: 0.500000 +DEAL::Integrating dealii::Functions::CutOffFunctionC1<2> -2,2 cube: 0.999764 +DEAL::Testing dim = 3 +DEAL::Center: 0.00000 0.00000 0.00000 +DEAL::Radius: 1.00000 +DEAL::Integrating dealii::Functions::CutOffFunctionC1<3> -2,2 cube: 0.999478 +DEAL::Center: 0.500000 0.500000 0.500000 +DEAL::Radius: 0.500000 +DEAL::Integrating dealii::Functions::CutOffFunctionC1<3> -2,2 cube: 1.00089 diff --git a/tests/base/function_cutoff_02.cc b/tests/base/function_cutoff_02.cc index 5540e56745..d2997c14d6 100644 --- a/tests/base/function_cutoff_02.cc +++ b/tests/base/function_cutoff_02.cc @@ -117,4 +117,8 @@ main() test<1, Functions::CutOffFunctionCinfty>(); test<2, Functions::CutOffFunctionCinfty>(); test<3, Functions::CutOffFunctionCinfty>(); + + test<1, Functions::CutOffFunctionC1>(); + test<2, Functions::CutOffFunctionC1>(); + test<3, Functions::CutOffFunctionC1>(); } diff --git a/tests/base/function_cutoff_02.output b/tests/base/function_cutoff_02.output index 2f3fe3da81..d694330994 100644 --- a/tests/base/function_cutoff_02.output +++ b/tests/base/function_cutoff_02.output @@ -62,3 +62,24 @@ DEAL::Integrating dealii::Functions::CutOffFunctionTensorProduct<3> -2,2 cube: 1 DEAL::Center: 0.500000 0.500000 0.500000 DEAL::Radius: 0.500000 DEAL::Integrating dealii::Functions::CutOffFunctionTensorProduct<3> -2,2 cube: 1.00779 +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.00009 -- 2.39.5