*
* 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.
static const unsigned int no_component = numbers::invalid_unsigned_int;
/**
- * Constructor. Arguments are the center of the ball and its radius.
- *
- * If an argument <tt>select</tt> 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<dim>::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<dim> = Point<dim>(),
+ const double radius = 1.,
+ 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,
/**
* Set the radius of the ball to @p r
- *
- * @deprecated Use set_radius() instead.
*/
virtual void
set_radius(const double r);
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;
};
+ /**
+ * 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<s$ is the distance to the center, and
+ * $s$ is the radius of the sphere. If vector valued, it can be restricted to
+ * a single component.
+ *
+ * @ingroup functions
+ * @author Luca Heltai, 2019
+ */
+ template <int dim>
+ class CutOffFunctionC1 : public CutOffFunctionBase<dim>
+ {
+ public:
+ /**
+ * Constructor.
+ */
+ CutOffFunctionC1(
+ const double radius = 1.,
+ const Point<dim> = Point<dim>(),
+ const unsigned int n_components = 1,
+ const unsigned int select = CutOffFunctionBase<dim>::no_component,
+ bool integrate_to_one = false);
+
+ /**
+ * Function value at one point.
+ */
+ virtual double
+ value(const Point<dim> &p, const unsigned int component = 0) const override;
+
+ /**
+ * Function values at multiple points.
+ */
+ virtual void
+ value_list(const std::vector<Point<dim>> &points,
+ std::vector<double> & values,
+ const unsigned int component = 0) const override;
+
+ /**
+ * Function values at multiple points.
+ */
+ virtual void
+ vector_value_list(const std::vector<Point<dim>> &points,
+ std::vector<Vector<double>> & values) const override;
+
+ /**
+ * Function gradient at one point.
+ */
+ virtual Tensor<1, dim>
+ gradient(const Point<dim> & 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 <tt>radius</tt> around
//////////////////////////////////////////////////////////////////////
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
}
+
+ template <int dim>
+ CutOffFunctionC1<dim>::CutOffFunctionC1(const double r,
+ const Point<dim> p,
+ const unsigned int n_components,
+ const unsigned int select,
+ bool integrate_to_one)
+ : CutOffFunctionBase<dim>(r,
+ p,
+ n_components,
+ select,
+ integrate_to_one,
+ integral_C1[dim - 1])
+ {}
+
+
+ template <int dim>
+ double
+ CutOffFunctionC1<dim>::value(const Point<dim> & p,
+ const unsigned int component) const
+ {
+ if (this->selected == CutOffFunctionBase<dim>::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 <int dim>
+ void
+ CutOffFunctionC1<dim>::value_list(const std::vector<Point<dim>> &points,
+ std::vector<double> & 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<dim>::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 <int dim>
+ void
+ CutOffFunctionC1<dim>::vector_value_list(
+ const std::vector<Point<dim>> &points,
+ std::vector<Vector<double>> & 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<dim>::no_component)
+ values[k] = val;
+ else
+ {
+ values[k] = 0;
+ values[k](this->selected) = val;
+ }
+ }
+ }
+
+
+
+ template <int dim>
+ Tensor<1, dim>
+ CutOffFunctionC1<dim>::gradient(const Point<dim> &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>;
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
test<1, Functions::CutOffFunctionCinfty>();
test<2, Functions::CutOffFunctionCinfty>();
test<3, Functions::CutOffFunctionCinfty>();
+
+ test<1, Functions::CutOffFunctionC1>();
+ test<2, Functions::CutOffFunctionC1>();
+ test<3, Functions::CutOffFunctionC1>();
}
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
test<1, Functions::CutOffFunctionCinfty>();
test<2, Functions::CutOffFunctionCinfty>();
test<3, Functions::CutOffFunctionCinfty>();
+
+ test<1, Functions::CutOffFunctionC1>();
+ test<2, Functions::CutOffFunctionC1>();
+ test<3, Functions::CutOffFunctionC1>();
}
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